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Kurzfassung 


Aufgrund ihres Leichtbaupotenzials bei relativ geringen Kosten gewinnen glas- 
faserverstärkte Polymere in industriellen Anwendungen zunehmend an Bedeu- 
tung. Sie verbinden die hohe Festigkeit von Glasfasern mit der Beständigkeit 
von z.B. duroplastischen Harzen. Bei der Verarbeitung von faserverstärkten 
Duroplasten kommt es zu einer chemischen Reaktion des Harzes. Die chemi- 
sche Reaktion geht mit einer chemischen Schrumpfung einher. In Verbindung 
mit der thermischen Ausdehnung kann das Material bereits beim Herstellungs- 
prozess beschädigt werden. Auch wenn das Komposit nicht vollständig ver- 
sagt, kann es zu Mikrorissbildung kommen. Diese Schäden können die Blast- 
barkeit des Bauteils und damit seine Lebensdauer beeinträchtigen. Faserver- 
stärkte Duroplaste enthalten Strukturen auf verschiedenen Längenskalen, die 
das Verhalten des Gesamtbauteils beeinflussen und daher für eine genaue Vor- 
hersage der Rissbildung berücksichtigt werden müssen. Das Verständnis der 
Mechanismen der Rissbildung auf den verschiedenen Längenskalen ist daher 
von großem Interesse. Auf der Grundlage von Molekulardynamiksimulationen 
wird ein Harzsystem zusammen mit einer Faseroberfläche und einer Schlichte 
auf der Nanoskala betrachtet und ein systematisches Verfahren für die Entwick- 
lung eines ausgehärteten Systems vorgestellt. Eine zweistufige Reaktion, eine 
Polyurethanreaktion und eine radikale Polymerisation, wird auf der Grundla- 
ge eines etablierten Ansatzes modelliert. Anhand des fertig ausgehärteten Sys- 
tems werden Auswertungen über gemittelte Größen und entlang der Norma- 
lenrichtung der Faseroberfläche durchgeführt, was eine räumliche Analyse der 
Faser-Schlichtharz-Grenzfläche erlaubt. Auf der Mikrolängenskala werden die 
einzelnen Fasern räumlich aufgelöst. Mit Hilfe der Kontinuumsmechanik und 
der Phasenfeldmethode wird das Versagen während des Aushärtungsprozesses 
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auf dieser Längenskala untersucht. In der Materialwissenschaft wird die Pha- 
senfeldmethode häufig zur Modellierung der Rissausbreitung verwendet. Sie 
ist in der Lage, das komplexe Bruchverhalten zu beschreiben und zeigt eine 
gute Übereinstimmung mit analytischen Lösungen. Dennoch sind die meis- 
ten Modelle auf homogene Systeme beschränkt, und nur wenige Ansätze für 
heterogene Systeme existieren. Es werden bestehende Modelle diskutiert und 
ein neues Modell für heterogene Systeme abgeleitet, das auf einem etablierten 
Phasenfeldansatz zur Rissausbreitung basiert. Das neue Modell mit mehreren 
Rissordnungsparametern ist in der Lage, quantitatives Risswachstum vorherzu- 
sagen, wo die etablierten Modelle eine analytische Lösung nicht reproduzieren 
können. Darüber hinaus wird ein verbessertes Homogenisierungsschema, das 
auf der mechanischen Sprungbedingung basiert, auf das neuartige Modell an- 
gewandt, was zu einer Verbesserung der Rissvorhersage selbst bei unterschied- 
lichen Steifigkeiten und Risswiderständen der betrachteten Materialien führt. 
Zudem wird zur Erzeugung digitaler Mikrostrukturen, die für Aushärtungssi- 
mulationen im Mikrobereich verwendet werden, ein Generator für gekrümm- 
te Faserstrukturen eingeführt. Anschließend wird die Verteilung mechanischer 
und thermischer Größen für verschiedene Abstraktionsebenen der realen Mi- 
krostruktur sowie für verschiedene Faservolumenanteile verglichen. Schließ- 
lich wird das neue Rissausbreitungsmodell mit dem Aushärtungsmodell kom- 
biniert, was die Vorhersage der Mikrorissbildung während des Aushärtungspro- 
zesses von glasfaserverstärktem UPPH-Harz ermöglicht. 
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Abstract 


Because of their lightweight potential at relatively low cost, glass-fiber rein- 
forced polymers are becoming increasingly important in industrial applications. 
They combine the high strength of glass fibers with the durability of, for exam- 
ple, thermoset resins. During processing of fiber-reinforced thermosets, the 
resin material undergoes a chemical reaction. The chemical reaction is ac- 
companied by chemical shrinkage. Combined with thermal expansion, the mi- 
crostructure can be damaged during this manufacturing process. Even if the 
component does not fail completely, micro-cracking can occur. This damage 
can affect the overall performance of the component and therefore its lifetime. 
Fiber-reinforced thermosets contain microstructures on different length scales, 
which influence the behavior of the overall composite and must therefore be 
taken into account for accurate prediction of crack formation. Understanding 
the mechanisms of crack initiation at different length scales is therefore of great 
importance. To this end, using molecular dynamics simulations, a resin system 
is considered together with a fiber surface and a sizing layer at the nanoscale. A 
systematic procedure for the development of a final cured system is presented. 
A two-stage reaction, a polyurethane reaction and a radical polymerization, is 
modeled based on an established molecular dynamics approach. For the fully 
cured system, evaluations are carried out in terms of averaged quantities as 
well as quantities along the normal direction of the fiber surface, resulting in 
spatially varying properties along the fiber-sizing-resin interface. At the micro 
scale, the individual fibers are resolved. Based on continuum mechanics and 
the phase-field method, the fracture during the curing process is studied in this 
work. In materials science, the phase-field method is widely used to model 
crack propagation. It is capable of describing complex fracture behavior and 
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shows good agreement with analytical solutions. Nevertheless, most models 
are introduced for homogeneous systems, and few approaches are available for 
heterogeneous systems. Existing models are discussed and a new model for 
heterogeneous systems is derived based on an established phase-field approach 
to crack propagation. The new multi-crack order parameter model is able to 
predict crack growth quantitatively, whereas the established models fail to re- 
produce an analytical solution. Furthermore, amore advanced homogenization 
scheme based on the mechanical jump condition is applied to the novel model, 
leading to an improvement in crack prediction even for a high contrast in stiff- 
ness and crack resistance of the considered materials. A curved fiber structure 
generator is introduced to generate digital microstructures that are used for cur- 
ing simulations at the micro scale. Subsequently, the distributions of mechan- 
ical and thermal quantities are compared for different levels of abstraction of 
the real microstructure, as well as for different fiber volume fractions. Finally 
the novel crack propagation model is combined with the curing model, which 
allows the prediction of micro-crack formation during the curing process of 
glass fiber reinforced UPPH resin. 
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1 Introduction 


1.1 Motivation and objectives 


Glass fiber-reinforced polymers are increasingly important in industrial appli- 
cations due to their lightweight potential at relatively low cost. They combine 
the high strength of glass fibers with the durability of, for example, thermosets. 
A novel unsaturated polyester polyurethane hybrid (UPPH) resin system al- 
lows co-molding of continuous fiber-reinforced polymers in e.g. compression 
molding. The relatively low cycle times allow a wide range of uses, such as 
in automotive applications [6, 7]. Understanding the behavior of such a com- 
posite, including manufacturing and possible damage during this process, is 
essential for designing safe and reliable products. 


Fiber-reinforced thermosets (FRTS) contain structures on of different length 
scales, cf. Figure 1.1. Starting with the nano length scale where individual 
atoms are resolved and the interfaces between fiber, resin and fiber sizing take 
place. A well-established method for this scale is molecular dynamics, for 
which a approach for investigating the fiber-sizing-resin interface during poly- 
merization will be presented in this work. At the micro length scale, the indi- 
vidual fibers are resolved. Based on continuum mechanics and the phase-field 
method, the fracture during the curing process will be studied. The meso length 
scale deals with whole fiber bundles, while the marco length scale represents 
the entire FRTS component. 
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Component 


Fibers Fiber bundles 


Fiber-sizing-resin 


length scale 


micro meso 
Phase-field method 


Figure 1.1: The different length scales of fiber-reinforced thermosets and their corresponding dif- 
ferent computational methods used in this work. With images reprinted with permis- 
sion from [8, 9]. 


Crack propagation modeling! Understanding failure and fracture behav- 
ior is a challenge in modern engineering and materials science, especially when 
considering the growing number of materials with a complex morphology and 
heterogeneous material properties, such as fiber-reinforced polymers (FRP). 
Predicting the resistance to failure and the complex crack propagation paths of 
such materials will improve the ability to determine effective load capacities 
and to develop efficient, safe, and predictable products. Linear Elastic Frac- 
ture Mechanics (LEFM) has proven to be capable of describing crack propa- 
gation in homogeneous materials in 2D [10, 11]. An extension to heteroge- 
neous materials is possible [12, 13], but a general approach which describes 
complex heterogeneous materials in 3D seems difficult and not feasible. An 
alternative approach is cohesive zone modeling (CZM), introduced by Baren- 
blatt [14] and Dugdale [15], which can be embedded in the finite element 
method (FEM), using cohesive finite elements. An overview of CZM can be 
found in Elices et al. [16]. Since these models demand conforming meshes, 


'The content of this section has been taken directly from Schiller et al. [1, 3] with minor linguistic 
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more advanced crack paths requires cumbersome remeshing methods. In com- 
parison, the generalized/extended finite element method (GFEM/XFEM) en- 
riches the solution space of the FEM to handle discontinuous functions [17]. 
This eliminates the need for conforming meshes and remeshing. However, both 
CZM and GFEM/XFEM are limited to describe complex crack propagation be- 
havior, including nucleation, branching, or the interaction between cracks. In 
a domain with sharp interfaces, different regions, e.g., phases or destroyed and 
unbroken material, occurring in the case of fracture, are distinctly separated, cf. 
e.g., Prahs and Böhlke [18], in the context of interface conditions on a sharp 
interface. This requires explicit tracking of the interface, which has proven 
impractical. 


Phase-field crack propagation models! An alternative approach to 
fracture utilizes the phase-field method (PFM), introducing order parameters 
to allow a smooth transition between various regions. This results in continu- 
ous order parameter fields, often referred as phasefields, and allows an implicit 
tracking of the interface on nonconforming meshes, and thus an efficient nu- 
merical treatment of singularities, such as grain boundaries or cracks. Thus, the 
PFM is widely established to describe the evolution of microstructures, such 
as solidification or solid-solid phase transitions, including consideration of dif- 
ferent types of physical aspects, e.g., thermodynamics, chemistry, or mechan- 
ics [19-23]. Phase-field approaches to brittle fracture have been developed in 
both the physical [24-26] and the mechanical community [27-29]. The latter 
is based on Griffith’s theory [30] and the variational formulation of Francfort 
and Marigo [31] and Francfort and Bourdin et al. [32]. Other more advanced 
applications, for example, deal with plasticity [33-36] or multiphysics [37-40]. 
For most of these models, the material is considered homogeneous. This is a 
reasonable assumption on macroscopic length scales. Often, however, failure 
mechanisms occur at smaller length scales, where many materials are hetero- 
geneous. Therefore, models that are able to describe fracture of heterogeneous 
systems are highly desirable. 
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Crack propagation in heterogeneous materials! Most phase-field 
models describing crack propagation in such systems introduce a varying crack 
surface energy. This is achieved either by anisotropy [41, 42] or interpolation 
of the surface energy [43, 44]. Schneider et al. [45] proposed a model that ex- 
tends the multiphase-field model of Nestler et al. [46], so as to describe crack 
propagation in polycrystalline systems. The model takes into account damage 
due to a phase transition to a common crack phase. This concept has also 
been extended to anisotropy [47] and plasticity [48]. A common approach to 
consider an interfacial fracture toughness is to lower the crack resistance in the 
interfacial region. Hansen-Dörr et al. [49, 50], for instance, model a locally 
varying value based on a virtual phase transition. To account for interfacial ef- 
fects, CZM have also been introduced into phase-field crack propagation mod- 
els [51-54]. Although these models can describe complex crack propagation 
in heterogeneous materials, including interfacial effects, and agree with the 
LEFM and experiments, they can lead to non-physical behavior as discussed 
by Henry [55]. 


An objective of this work is to introduce a novel multi-crack order parameter 
(MCOP) phase-field model for fracture which is able to overcome issues of 
the established single-crack order parameter (SCOP) approaches for modeling 
crack propagation in heterogeneous systems based on the phase-field method. 
Therefore, some issues of the SCOP model are illustrated in simple 1D and 
2D simulation setups and several advantages of the novel model are demon- 
strated. Moreover, the work at hand is to extend this model by incorporating 
a more sophisticated scheme for the underlying homogenization problem [56- 
61]. Therefore, the work of Schneider et al. [60] is applied to the MCOP model. 
Finally, an exemplary FRP system is used to demonstrate the limitations of the 
basic scheme and the advantages of the proposed model, where the homoge- 
nization scheme is based on mechanical jump conditions. 


1.1 Motivation and objectives 


The role of fiber sizing? Besides the matrix material and the fiber, the siz- 
ing (fiber coating) plays a crucial role in the manufacturing and performance 
of FRPs. The fiber size consist of multiple components and fulfills a variety 
of tasks, such as protection of the fiber, improvement of the fiber handling, or 
enhancement of the adhesive bonding of fiber and matrix [62]. Early investi- 
gations considering the sizing were conducted by Plueddemann [63, 64] and 
Loewenstein [65]. They mainly focus on the coupling agent, which provides 
the functionality of bonding to both the glass fiber and the resin. In the litera- 
ture, Thomason [66-68], Gao and Mäder [69] and Liu et al. [70], among others, 
examined the sizing of glass fibers with different sizing formulations and using 
various approaches. Furthermore, Karger-Kocsi et al. [71] summarize recent 
advances in interphase technology for several fibers, and matrix materials, as 
well as sizes. Also, Thomason [62] recently gave a detailed overview of glass 
fiber sizing. He pointed out that the actual sizing is always a proprietary secret 
that leads to a lack of understanding of the fiber sizing, especially since the 
knowledge is very fragmented. As the fiber interphase is a common point of 
failure in FRPs, through mechanisms such as fiber pull-out and fiber debond- 
ing, this lack of knowledge is a severe impediment to improving these materi- 
als. Moreover, fiber interphase has been mostly studied in the literature based 
on an experimental approach. In contrast, few investigations have been con- 
ducted based on simulative approaches. Therefore, modeling of the fiber-resin 
interphase, including the sizing, e.g. based on molecular dynamics, is highly 
desirable as it provides insights into the interphase and could improve the un- 
derstanding of FRP. In particular, regarding their manufacturing, performance, 
or failure mechanism, this could lead to an improvement of this whole class of 
materials. 
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Only a few attempts to investigate the interface between fiber and matrix, based 
on molecular dynamics (MD), can be found in the literature. While the inter- 
phase between carbon fibers was examined [72-74], including the fiber siz- 
ing [75], no investigation of the interphase for glass fibers and their sizing have 
been conducted so far. Insight into the formation of the network during the 
polymerization, in combination with the fiber size, could improve the under- 
standing of sizing extensively. 


Curing of FRTSs During the processing of FRTS, e.g. compression mold- 
ing, the resin material undergoes a chemical reaction. The formation of the 
final polymer structure gives the component its final strength and properties, 
such as mechanical stiffness. Accompanying the chemical reaction, the curing 
process, is a chemical shrinkage. Together with a thermal expansion due to 
the higher temperature, the microstructure can already exhibit stresses during 
this manufacturing process. In particular, the high contrast between the ther- 
mal, chemical and mechanical material behavior of fiber and matrix can cause 
damage at this stage. Even if a complete failure of the component does not 
occur, microcracks may occur, e.g. at the fiber-resin interface. This damage 
can affect the overall performance of the component and therefore its service 
life. Understanding the mechanism of crack formation at the length scale of 
individual fibers is therefore of great interest. 


In this work, a structure generator for curved fibers is introduced. Using such 
volume elements and material properties derived from the study of the resin 
system with molecular dynamics simulations [76], curing simulations are per- 
formed. Subsequently, the distribution of mechanical and thermal quantities 
will be compared for different levels of abstraction of the real microstruc- 
ture as well as for different fiber volume fractions. Finally, the combination 
of the novel MCOP model with the mechanical jump condition framework, 
cf. eg. [60], allows the prediction of microcrack formation during the curing 
process of glass fiber-reinforced UPPH resin. 


1.2 Notation and conventions 


Objectives The main objectives and novelties of this work are listed in the 
following: 


1.2 


The two-stage polymerization process of the UPPH resin is studied by 
molecular dynamics simulations. In contrast to Schwab and Dennis- 
ton [76], the glass fiber and sizing are also considered, and the evolution 
of various quantities during polymerization along the fiber-sizing-resin 
interface is examined. 


A volume element generator for long curved fibers is developed based 
on molecular dynamics and the discrete element method. 


The curing process of UPPH in such glass fiber-reinforced volume el- 
ements is studied and stresses due to thermal and chemical strains are 
discussed. 


A novel phase-field crack propagation model based on multiple crack 
order parameters is developed. It is demonstrated that it provides qual- 
itatively and quantitatively better predictions of crack formation in het- 
erogeneous materials than established models. 


In the novel MCOP, the mechanical jump condition framework [22] is 
incorporated. The advantages of such an extension are shown, e.g., by 
crack formation around a single inclusion problem. 


Crack propagation in FRTS is studied. External displacement boundary 
conditions are used to demonstrate the advantages of the MCOP as well 
as the jump condition framework. Finally, crack initiation during curing 
of glass fiber-reinforced UPPH resin is predicted using the novel model. 


Notation and conventions 


In this work, a symbolic notation for tensors is used: Scalar-valued quantities 


are denoted by a normal font, e.g. a,b. Furthermore are 1" and 2" order 
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tensors, or vectors, denoted by bold symbols, e.g. a, B. In addition, 4" order 
tensors are denoted by double-strike symbols, e.g. C. The tensors are restricted 
to an orthonormal basis (ONB) {ex, ey, ez} with the unit vectors e; Vi = x,y, z. 
Greek letters, e.g. &,ß, are used to denote phase indices in the context of the 
phase-field model. 


Operations The scalar product is denoted by a -b and A- B for vectors a, b 
or 2" order tensors A, B, respectively. A linear mapping of a 2" order tensor 
A and a vector bis defined by Ab, and for a 2™ order tensor A and a 4" order 
tensor C by C[A]. The dyadic product of two vectors follows by a® b. The 
vector norm can be written with |a| = v/a - a, and the tensorial Frobenius norm 
analogously with |A| = VA- A. The inverse of a tensor A is denoted by A`! 
and the transposed tensor by AT. The trace of a 2™ order tensor A is defined 
by tr(A) = A-1, where the 2™ order identity tensor is 1. 


Derivatives Total derivatives are denoted by de a for a vector a with respect 
to x, while partial derivatives are denoted by 2a, The derivatives respectively 
apply to any tensorial order of the numerator or denominator. For the tempo- 
ral derivative the abbreviation @ = ap for an arbitrary field @ and time t is 
introduced. The spatial gradient æ with respect to a position vector x can 
be written by grad(@) or using the nabla operator with Va. The divergence 
of a vector field a is denoted by div(a) = grad(a) -1 or V -a and can also 


be extended to higher tensorial orders. The Laplacian operator is defined by 
Ag = div(grad(9)). 


1.3 Outline of the thesis 


1.3 Outline of the thesis 


This dissertation is organized as follows: Selected fundamentals of compos- 
ite materials, continuum mechanics, phase-field modeling, fracture mechanics 
and molecular dynamics are introduced in Chapter 2. Subsequently, in Chap- 
ter 3, the two-stage polymerization of the fiber resin system considered is in- 
vestigated using molecular dynamics. In order to study this curing process 
at the micro scale, a structure generator for curved fibers in FRPs is estab- 
lished in Chapter 4. Based on the generated structures, micro-scale curing 
simulations are performed and the corresponding curing model is formulated 
in Chapter 5. In order to be able to describe crack propagation in such a het- 
erogeneous system, different crack propagation models based on phase-field 
modeling are presented in Chapter 6. These are qualitatively and quantitatively 
compared, resulting in a novel crack propagation model for heterogeneous sys- 
tems, including advanced homogenization schemes. The ability of the model 
to predict crack propagation in FRPs is demonstrated and finally applied to the 
compression molding process of glass fiber-reinforced UPPH resin in Chap- 
ter 7. Thereby, the loading due to thermal expansion and chemical shrinkage 
on the micro-length scale is considered. Finally, conclusions and a outlook are 
given in Chapter 8. 


2 Selected fundamentals 


2.1 Composite materials and compression 
molding process 


A composite material consists of two or more components. Typically, the com- 
ponents have significantly different material properties, such as chemical or 
mechanical properties and are often composed of materials of different classes, 
such as steel-reinforced concrete, fiber-reinforced ceramics, or fiber-reinforced 
polymers. The general idea of a composite material is to combine the con- 
stituents to get a composite material with better effective properties. For exam- 
ple, for fiber-reinforced polymers (FRP), the matrix materials primarily embed 
the fibers, transfer loads and protect the fibers. In contrast, the fibers mainly 
provide the strength and stiffness of the composite. The effective strength is 
significantly increased compared to the pure matrix materials, while other prop- 
erties, e.g. chemical resistance, may be reduced compared to a pure glass. 


In this work, FRPs, in particular fiber-reinforced thermosets (FRTS), are con- 
sidered in the context of the compression molding process, which is briefly 
introduced in the following. The reader is referred to [7, 77] for a more com- 
prehensive introduction to FRTS and the compression molding process. 
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Resin 


Chopped 
fibers 


Upper carrier foil 


Impregnation 


Fiber rovings 


Lower carrier foil 


Semi-finished SMC 


Figure 2.1: Manufacturing process of semi-finished sheet molding compound (SMC). Reprinted 
from [78]. 


Glass fiber-reinforced UPPH!' Fiber-reinforced polymers are increas- 
ingly important in industrial applications, including aerospace, automotive, 
marine, and construction industries. The majority of these are glass fiber- 
reinforced polymers, probably due to the relatively low cost while still offering 
good performance. A relatively new thermoset matrix material in this con- 
text is the unsaturated polyester polyurethane hybrid (UPPH) resin system. It 
combines a polyurethane (PU) polyaddition with a radical polymerization of 
unsaturated polyester (UP), resulting in a two-step reaction of the thermoset. 
The hybrid networks formed by this copolymerization increase overall prop- 
erties such as crack resistance, tensile strength and toughness [6]. There is 
therefore increasing interest in using this material in industrial applications 
such as the automotive industry [7]. 


Compression molding process For FRPs, the choice of manufacturing 
process is often based on the fiber aspect ratio, i.e., the ratio of fiber length 
to diameter. Depending on the ratio, different high shear stresses can occur 


'The content of this section has been taken directly from Schiller et al. [2] with minor linguistic 
changes. 
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i : Heed tt TEREST 
A ER ; : 
Initial charge positioning Mold flow Resin curing 


Figure 2.2: The compression molding process: Placing the SMC initial charge, compression mold- 
ing and finally, curing of the resin. Reprinted from [78]. 


during the processes, which prevents a free choice of the manufacturing pro- 
cess. In this work, the compression molding process is chosen because it 
also allows the co-molding of continuous fiber-reinforced polymers as local 
reinforcements, leading to continuous-discontinuous FRPs. This class of com- 
posites combines the possibility of complex structures of discontinuous fiber- 
reinforced polymers with the high strength and stiffness of continuous fiber- 
reinforced polymers and offers an immense weight saving potential. In this 
context, compression molding allows low cycle time and therefore low cost. 


The compression molding process of FRTS uses a prepreg, a sheet molding 
compound (SMC). To produce the SMC, the resin is mixed and applied evenly 
to two films through doctor boxes. Chopped glass fibers are then dropped 
randomly onto the films. The films are then joined and the resin cures during 
a polyurethane reaction to a so-called B-stage, cf. Figure 2.1. The finished 
SMC can be easily transported and cut to size. After the films are removed, the 
charge is placed in the press onto the hot molds. By closing the molds, pressure 
and heat are applied and the SMC charge flows into its final shape and radical 
polymerization takes place, curing the part to its final strength, cf. Figure 2.2. 
After the molds are opened, the part is removed and cooled to room temperature 
for further processing. An example of the final microstructure is shown in 
Figure 2.3. The chopped fiber roving, consisting of hundreds of individual 
fibers, can be identified in the CT scan. Based on the flow direction, the bundles 
may have a preferred direction. Thus, a planar distribution is observed in the 
CT scan due to the planar geometry of the final component. 
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Figure 2.3: 3D visualization of the tracked and clustered SMC fiber bundles. Only clustered fiber 
bundles with a minimum of 100 sub-bundles are visualized. Reprinted with permission 
from [9]. Copyright 2021 Elsevier. 


2.2 Continuum mechanics 


This section gives a brief introduction to continuum mechanics, mainly fol- 
lowing Altenbach [79] and Bertram and Gliige [80]. For more detailed funda- 
mentals of of continuum mechanics, the reader is referred to e.g. [79-84]. In 
continuum mechanics, material bodies are considered which continuously oc- 
cupy a region in three-dimensional Euclidean space R? by an infinite number 
of material points. This work is restricted to Cauchy-Boltzmann continua, for 
which each material point has three translational degrees of freedom and, for 
example, no rotational degrees of freedom. 


2.2.1 Kinematics 


A material body Q with a boundary ð Q and a unit normal vector n directed 
outwards is considered. Any material point within Q can be described by a 
position vector X at a reference time tọ in a reference configuration. The 
position vector at a time f > tọ can be described by 


x =X(X,t) (2.1) 
in the current configuration by the motion function 
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Deformation The gradient of the motion x with respect to the reference 
configuration 


2 (2.2) 


is called the deformation gradient and describes the transformation of a line 
element dX of the reference configuration into a line element dæ in the current 
configuration by 

de=FdX. (2.3) 


In addition, the transformation of an area element dA and a volume element 
dV of the reference configuration into elements of the current configuration 
(da,dv) follows by 


da = det(F) (F') dA, (2.4) 


dv = |det(F)| dV. (2.5) 


The difference between the reference and current position vector yields the 
displacement vector 
u=X(X,1)-X, (2.6) 


and time derivative of the motion x the velocity vector 


mS Ox(X,t) 


oe (2.7) 


The gradient of the displacement vector with respect to the reference configu- 
ration 


= —_=F-1 (2.8) 


is called the displacement gradient and can also be directly related to the defor- 
mation gradient with the identity tensor 1. 
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Strain measures Based on a polar decomposition of the deformation gra- 
dient the symmetric Green’s strain tensor 


ES = ; (F'F-1) (2.9) 


can be defined, which is independent of of rigid body motion, and therefore is 

a common strain measure for finite deformation. The symmetric infinitesimal 
strain tensor 1 

=-(H+H" 2.10 

=, (H+H"), (2.10) 


is the linearized variant of ES, which is only independent of translational rigid 
body motion, not to rotational rigid body motion and is therefore often used in 
the context of small deformations. 


2.2.2 Balance laws 


This section considers sufficiently smooth physical fields without singular sur- 
faces. For a discussion of continua with material singular surfaces, the reader 
is referred to e.g. Prahs and Böhlke [18]. 


Balance of linear momentum The balance of linear momentum can be 


written by 
d 
F) podv= | pbav+ | tda, (2.11) 
dt Jo Q dQ 


with the mass density p, and implies that the change of linear momentum 
Jo pvdv is equal to the forces of a body force vector b and the traction vec- 
tor t. 


The Theorem of Cauchy states that the traction vector is given by 


t=on, (2.12) 
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with the normal vector n and the Cauchy stress tensor o [80]. For the static 
case and when no body force is present, the local balance of linear momentum 
can be written by 

div(o) =0. (2.13) 


Furthermore, the balance of angular momentum is satisfied with the symmetry 
of the stress tensor 
o=a', (2.14) 


as stated by Boltzmann’s axiom [80]. 


First law of thermodynamics The first law of thermodynamics, or the 
balance of total energy, can be defined by 


| Le(et zee) av= [portna+ | (t-v—q-n)da, (2.15) 
dt Ja 2 Q dQ 


with the specific internal energy e, a specific heat source r and the heat flux 
vector q. Therefore, the change in total energy, the sum of internal and kinetic 
energy, is equal to the power of a volumetric and surface forces and the heat 
supplied by a heat source and flux. 


Second law of thermodynamics The first law of thermodynamics pro- 
vides a balance of total energy, but it does not state how energy can be con- 
verted, nor in what direction. Therefore, the second law of thermodynamics is 
given by 


d r q 
= d dv+ -nda > 2.1 
le m 0 


with the absolute temperature 0, the specific entropy 17, and states that the en- 
tropy production of a body is never greater than its entropy change due to the 
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heat flux and the heat source. In combination with the first law of thermody- 
namics, the localized Clausius-Duhem inequality (CDI) 


o-é—pw—pnd >20 (2.17) 


with the Helmholtz free energy y = e — On and the temperature gradient g = 
00/0X can be formulated. This inequality equation imposes restrictions on 
the constitutive equations describing the material response, e.g. by following 
the work of Coleman and Noll [85]. 


2.2.3 Infinitesimal deformation 


The deformation of a body is infinitesimal, or small, if 
|\H|</|1| (2.18) 


is fulfilled. And therefore, e.g. the strain tensor can be approximated by the 
geometrically linearized variant 


ES xe. (2.19) 


Furthermore, for the gradient of an arbitrary field y with respect to the refer- 


ence and current configuration, 


oy ow 


IN Jr P ~ Po (2.20) 


applies. 
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2.2.4 Material theory 


The balance of linear momentum (2.13) contains nine unknown physical fields, 
three displacement and six stress components, but only provides three equa- 
tions. Therefore, six additional constitutive equations have to be formulated to 
relate the material response, hence the stresses, to the kinematic measures. In 
addition, some basic principles about the material behavior are postulated in 
the following [80]. 


e Principle of Determinism: The stresses at a material point are defined 
only by the displacement field of the body in the past and in the present. 


e Principle of Invariance under Superimposed Rigid Body Motions: 
The stresses are not directly caused by translations or rotations of the 
body. 


e Principle of Local Action: The stresses in a material point are deter- 
mined by the motion of only a finite neighborhood of that point. 


Hooke’s law The materials in this work are assumed to be elastic, which 
describe materials for which the current stresses depend only on the current 
deformations [80]. In the context of small deformations, the linear mapping of 


strain to stress tensor 
o =Clel, (2.21) 


with the 4" order stiffness tensor C, is the Hooke’s law in a general anisotropic 
form. The strain energy density follows by the quadratic form 


fa = ic [e]-e, (2:22) 


and the stress and stiffness can also be described by the potential relations 


= d fel = 3? fa = do 
en Tr ae (2.23) 


19 


2 Selected fundamentals 


In the isotropic case, Hooke’s law can be simplified to 
o =Atr(e)1+2pe, (2.24) 


with the Lamé parameters A and u, but can also be related to other elastic 
moduli, e.g. the Young’s modulus E and the Poisson’s ratio v, cf. e.g. [80]. 


Spectral form and principal stresses The symmetric stress tensor can 
be rewritten in the spectral form 


3 
o= y o'p'®p', (2.25) 


with the principal stresses oÝ Vi = 1,2,3 and the corresponding principal stress 
axes p’. The numbering of the principal stresses usually follows the convention 

0.080. (2.26) 
Stress deviator and von Mises stress In addition to the spectral form, 


the stress tensor can also be decomposed into a spherical part o° and the stress 
deviator o’ through 


o’=<-tr(o)l, oe =0-0°. (2.27) 


Note that for a hydrostatic stress state o = — pl at a pressure p, the deviator part 
vanishes, i.e., o’ = 0. Based on the stress deviator, the von Mises equivalent 


3 
OM = je lo’) (2.28) 


can be defined, which relates the three-dimensional stress state to a fictive uni- 


stress 


axial stress state. 
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2.3 Phase-field modeling 


The phase-field model (PFM) is a well-established approach for describing mi- 
crostructure evolution on the micro and meso length scales. Ina PFM a singular 
interface is described by a diffuse smooth transition instead of a sharp interface. 
This avoids a cumbersome explicit tracking of interfaces, which has proven to 
be unfeasible. In recent years, PFM has been widely used to describe the evo- 
lution of microstructures, such as solidification or solid-solid phase transitions, 
taking into account different types of physical aspects, e.g. thermodynamics, 
chemistry or mechanics [20-23]. For a comprehensive overview, the reader is 
referred to e.g. [19]. 


In this section, a basic multiphase-field model is presented based on the work 
of Steinbach et al. [86] and Nestler et al. [46]. Many of the extensions to the 
models that are state of the art nowadays are omitted, since the focus of this 
work is primarily on the phase-field model in the context of crack propagation. 
Instead, the model is solely used to parameterize the underlying morphology 
in a diffuse context. For an introduction to the PFM for crack propagation, the 
reader is referred to Chapter 6. 


2.3.1 Free energy 


For a material body © with a singular interface S, e.g. between two subregions 
with different material properties, a free energy functional 


F [u] = f Aula) dv-+ | frda (2.29) 
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is introduced. Thus, fpux describes the energy densities of the bulk regions 
which depends e.g. on the displacement vector u. The interfacial energy den- 
sity fintf, €-g., consists of the surface energy density of the surface. With the 
N-tuple & of order parameters and their corresponding gradients 


p= {o',0°,...,0%}, Ve = {Vo',V¢’,..., Vor}, (2.30) 


the free energy can be parameterized by these order parameters, yielding 


Flp, Vp, u] = hi (Fou (,V,t4) + fine, VB) dv = Í fi. (2.31) 


When the order parameters & are associated with an indicator functions, both 
functionals can be considered equivalent. In contrast, for a diffuse and smooth 
interface, where the order parameter can be considered as the local volume 
fraction of each subregion, the free energy (2.31) is an approximation of (2.29). 
The interfacial energy density can be further decomposed by 


Fintf Zu Fpot + Serad (2.32) 


into a potential contribution fpot and a gradient contribution fgrad. For the latter 
various formulation can be found in literature. In the work at hand 


faa (VO) = N) T Yap VO- VOP, (2.33) 


a B>a 


from Steinbach and Pezzolla [87] is used. Here, €¢, is a measure of the width of 
the diffuse interface and Ygg is the surface energy density between a subregion 
of a and B. For fpot, common choices in the literature are multi-well or multi- 
obstacle potentials. In this work, the latter is used in the form of 


foot( am 2 > Yap O° 9? (2.34) 
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from Nestler et al. [46]. In addition, the Gibbs simplex 


a-fo 


must be satisfied. The bulk energy density fouk in this work is limited to the 


N 
Der 1g z0va =i. (2.35) 


effective strain energy density, which can be decomposed into phase-specific 
strain densities {I by 


N 
foulo Vo, u) = ff =} AOSS, (2.36) 


with an interpolation function h*(&), which can be chosen arbitrarily with 
respect to some restrictions. The most straightforward choice is the order pa- 
rameter itself, h” (ġ) = @%, which will be used in the following work. 


2.3.2 Evolution equation 


To obtain an evolution equation for the order parameters @, a Ginzburg-Landau 
type evolution equation [88] is introduced. In the context of PFM it is often 
referred to as the Allen-Cahn evolution equation, cf. Allen and Cahn [89]. 
Nestler et al. [46] introduced in addition a Lagrange parameter to account for 
the sum constraint of (2.35), which yields 


ov = 59a tN 3 508 (2.37) 


M| F 1,8% 
Ep, i 


for the temporal evolution of the order parameter of the subregion &, with a 
mobility M . Steinbach and Pezzolla [87] introduced an alternative formulation 


; 1 Ww ÔF SF 
a YỌ mh u 3: 
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where dual interactions are considered, allowing for individual mobility M ap 
for each interaction. 


2.3.3 Considering mechanical jump conditions? 


The phase-specific strain energy densities are modeled by hyperelastic po- 


a 
tentials. Assuming small deformations 


1 
fi = a ‚et, o” =C* fe], (2.39) 


follows, with the phase-specific infinitesimal strain tensor e and the stiffness 
tensor C“. In general, the relation between phase-specific strains and the total 
strain tensor 


g= L o%e% = sym (grad (u)), (2.40) 


with the displacement vector u, is unknown. Using a Voigt-Taylor homoge- 
nization scheme, €% = e Va = 1,...,N is assumed. This results in a linear in- 
terpolation of the phase-specific parameters, e.g. stiffness. In contrast, Schnei- 
der et al. [60] presented a scheme which accounts for the mechanical jump 
conditions 


[H] aB aP g yo le] «B 08 = o, (2.41) 
with [y] = we yê, 


which was recently also applied to single-crack order parameter phase-field 
models by Prajapati et al. [47] and Hansen-Dörr et al. [90]. Thereby, the first 
equation resembles a kinematic compatibility, since the deformation gradient 
H can only exhibit a jump a®P in normal direction of the singular surface. In 
contrast, the balance of linear momentum on a material singular surface, c.f., 


>This section is based on the work of Schiller et al. [3]. Minor linguistic changes and additions 
have been made. 
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Figure 2.4: Schematic representation of a diffuse interface when mechanical jump conditions are 
considered. Each point of the interface can be associated with a laminate-like structure 
of different volume fractions of both phases and a normal direction n. 


e.g., Prahs and Böhlke [18], given by Equation (2.41)b, prohibits a jump of 
the stress vector normal to the singular surface. Regarding a multiphase-field 
approach, the normal vector n°8 of the singular surface between phase œ and 
B and the jump of the infinitesimal strain tensor, regarding the diffuse interface 
context is given by 


aß vo -Vof 


ze 2), 2.42 
|vo@—voß| 2 


le]? = sym (aP Q n8) : 


In order to solve the governing equations, the unknown jump vectors aß have 
be to determined. Reformulating the problem as a system of linear equations 
allows us to determine an effective stiffness for fixed order parameters [60, 90]. 
Based on a staggered approach, the governing equations are solved iteratively 
with an additional static criterion for crack propagation [1]. A more detailed 
introduction to mechanical jump conditions, in the context of a multiphase- 
field approach, is given by Schneider et al. [60]. 
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The introduction of order parameters and the interpolation of phase-specific pa- 
rameters requires additional assumptions. This is extensively discussed in the 
literature, e.g. by Ammar et al. [58], Durga et al. [56], or Mosler et al. [57], 
among others. Basic schemes such as a linear or harmonic interpolation of the 
phase-specific stiffnesses represent the upper and lower bounds of the strain 
energy in the context of a homogenization theory. Schneider et al. [22] dis- 
cussed the ability to reproduce sharp interface properties by a multiphase-field 
model when considering mechanical jumps. Figure 2.4 depicts a schematic 
representation of such diffuse interface when mechanical jumps are consid- 
ered. Therefore, each point of the interface can be thought of as a laminar 
structure. Instead of only the local volume fractions, also the normal vector n 
of the interface determines the effective elastic material behavior. In addition, 
a crack propagation model with multiple crack order parameters is introduced 
in this work, cf. Chapter 6. For such a model, the ability to obtain an accurate 
approximation of the phase-specific energy densities for sharp interfaces is of 
even higher interest. 


2.4 Linear elastic fracture mechanics 


This work is limited to crack propagation in the context of linear elastic frac- 
ture mechanics (LEFM), i.e., the material exhibits only elastic behavior. In- 
stead, the inelastic processes that may be associated with crack growth, such 
as microstructural defects, plasticity, or void formation, are considered within 
the process zone. This zone, around the crack tip, is assumed to be small com- 
pared to the length scale of the body. Note that this zone is not resolved in the 
context of LEFM and this work. This section contains only the most necessary 
basics of LEFM, for a more comprehensive discussion the reader is referred to 
e.g. [91]. 
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I II II 


Figure 2.5: Different crack modes: Mode I has the crack opening normal to the crack plane. And 
Mode II and Mode III with the crack opening tangential to the crack plane in different 
directions. 


Crack opening modes When considering a crack with a crack tip and a 
crack surface, three basic crack opening modes can be identified. For mode I, 
the crack opening is normal to the crack surface and can be associated with a 
tensile stress normal to the crack plane. For both Mode II and Mode III, the 
crack opening is tangential to the crack plane. Mode II can be associated with 
a sliding mode and with shear stress parallel to the crack plane and perpendic- 
ular to the crack front. Mode III can be considered as a tearing mode, with 
shear stress parallel to the crack plane and front. These three modes are also 
displayed in Figure 2.5. Note that for most real microstructures and loads a 
mixed mode, a superposition of the three modes, can be expected. 


Stress intensity factors For pure crack opening modes, the fields near the 
crack tip, e.g. stresses and strains, can be fully characterized by stress intensity 
factors Kj, Kj or Ky for each corresponding crack mode. These factors depend 
on the geometry of the body, the crack geometry and the loading and can be 
determined analytically for some cases, cf. eg. [91]. Based on the K-concept, 
which assumes that the crack tip and the process zone can be described by the 
stress intensity factors, the criterion 


Ki = Kic, Ky = Kite; Km = Kine; (2.43) 
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for crack growth of the corresponding pure crack mode can be formulated. The 
fracture toughness values Kj, Vi = I, II, II are material-specific properties. In 
the case of a mixed mode, an equivalent criterion f(K1,Ku,Kin) = 0 can be 
formulated with a generalized function f [91]. 


Energy release rate Based on the energy potential of a material body IT, 
including the potential of external forces, the energy release rate G follows 
from 
dII 
G=- aa (2.44) 
where dA is an infinitesimal crack surface. With the energy release rate, the 
criterion 


G=G, (2.45) 


can be formulated for crack growth with crack resistance or critical energy 
release rate G.. A more general approach to this criterion is Griffith’s fracture 
criterion [30] 


—-+—=0, (2.46) 


with the fracture energy I‘. From Equation (2.44) it follows that the energy 
released by crack growth must be equal to the energy required to create the cor- 
responding new crack surface, so that the crack resistance can also be related 
to the surface energy density y of a crack flank by G. = 2y. 


2.5 Molecular dynamics 


General idea? Molecular dynamics (MD) is a simulation method that al- 
lows prediction of the temporal evolution and interaction of atoms and mol- 
ecules, mostly using numerical integration of Newtons equations. The MD 


3The content of this section has been taken directly from Schiller et al. [2] with minor linguistic 
changes. 
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method has its origins in the work of Alder and Wainwright [92, 93] and Still- 
inger and Rahman [94]. An overview of the history of MD and its application 
is given by van Gunsteren and Berendsen [95]. 


2.5.1 Lagrange equations 


The structure of the derivation of the equations for MD is based on the La- 
grangian formalism, cf. Goldstein [96]. Based on the action functional 


S= | Lan... ano... ow)d (2.47) 


with the position vector æ; for a particle i and the velocity vector v; = &; for N 
particles in the system considered. The Lagrangian follows by 


L (£1, .., EN, U1,- UN) =K(vı,...,un)-Uleı,...,En); (2.48) 


with the kinetic energy K and the potential energy function U. For a stationary 
point of S, 


d (IL) IF 
era 2 


must hold, which are the Euler-Lagrange equations in the context of molecular 
dynamics. The kinetic energy K is given by 


1 N 
K(v1,... 0N) = 5 Lmivi- vi, (2.50) 


with the mass m; of particle i. The potential energy function U can be associ- 
ated with the conservative force vector f; for an particle i with 


d 
Fila, ‚an)= =U (21,---, 2N). (2.51) 


Ti 
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Algorithm 1: Velocity Verlet algorithm 


t< to 
while t < feng do 
1 
x; (t+At) xi (t) + Aty; (t) 4 „Ara; (t) 
1 1 
Vi (+ zar) -v(t)+ „Ara; (t) 
i (t+ At 
aj(t+At)h + getty Bt) 
m; 

1 1 
vi(tt+At) Cv; t+ 54t + 5Ata; (t+ At) 
t<t+At 

end 


With (2.50) and (2.51) the Equations (2.49) can be rewritten by 
ma; — fi=0, (2.52) 


with the acceleration vector a; = 0; = &;, which resembles Newton’s second 
law of motion and is the basis of the molecular dynamic simulations performed 
in this work. 


2.5.2 Microcanonical ensemble 


A system with N particles, e.g. atoms and additionally a fixed volume of the 
system (V) and a fixed internal energy (E) is considered as a microcanonical 
ensemble. This is often referred to as a (constant) NVE ensemble. Applying 
a time integration method to the equations of motion (2.52) yields a molecular 
dynamics algorithm. 
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Verlet integration Based on Newton’s equations of motion, a Taylor ex- 
pansion of the position coordinates x; can be performed at two different time 
steps, t + At and t — At. Combining these expansions yields 


x; (t+At) = 2a; (t) — x(t — At) + l Ara; ()+0(Ar*), (2.53) 


with the time step width At and an error term of order Ar*. This integration 
scheme is also referred to as Verlet integration [97]. 


Velocity Verlet integration The Verlet integration evolves solely the posi- 
tion vector. Although the velocity could be constructed at any time, an alterna- 
tive integration scheme is given by 


1 
xi (t+ At) = x; (t)+Aty; (t+ At) + „Ara; (t+At) (2.54) 


a (1) + 5 At (a; (eee ean). (2.55) 


which is known as the Velocity Verlet integration [98] and evolves the position 
as well as the velocity vector. Despite the rather simple scheme, this integrator 
exhibits some important properties, such as time reversibility, and is therefore 
used as the time integration scheme for the molecular dynamics simulations of 
this work. A possible implementation of this scheme is shown in Algorithm 1. 


2.5.3 Isobaric-isothermal ensemble 


As an alternative to the microcanonical ensemble, also the number of particles 
(N), the pressure (P) and temperature (T) can be fixed, resulting in the isobaric- 
isothermal or NPT ensemble. One way to archive this ensemble is to use a 
Nose-Hoover thermostat. Instead, in this work a NVE ensemble is used with 
a velocity Verlet integration, along with a thermostat and a barostat. For the 
latter, a Langevin thermostat [99] and a Berendsen barostat [100] are chosen, 
which are briefly introduced below. 
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Langevin thermostat The basic idea of the Langevin thermostat [99] is to 
add a random force ff to the particles with 


r kgm;0 
RS ay 2.56) 


the Boltzmann constant kg, and a damping factor y. This force mimics the 


force caused by random collisions of solvent atoms with the particles at a given 
temperature 0. In addition, a friction force 
f mi 


f= si (2.57) 


is added to damp the velocity of the particles. With forces of this nature it can 
be shown that the local temperature evolves towards 0. 


Berendsen barostat The Berendsen barostat [100] couples the system 
weakly to a heat bath of some pressure po. This adds an additional term 


dp| _Po=Pp 
dt | path Tp 


(2.58) 


describing the coupling to the heat bath, with the pressure of the system p and 
a time constant T,. This can be archived by rescaling the simulation domain by 


the factor 
KrAt 


Tp 


u=1- = (p-p), (2.59) 


at each time step with an isothermal compressibility Kr. 


2.5.4 Potential energy function 


The potential energy function U of a system with N particles can be decom- 
posed into 
U= pinta a unter (2.60) 
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where U' represents intramolecular interactions and U!" represents inter- 
molecular interactions. 


Intramolecular interactions The intramolecular part of potential energy 
function U'™"* can be further decomposed with 


yita BER L u) + us (2.61) 


i>j 


in a pairwise interaction contribution u) and an additional Coulomb potential 
US. Note that the subscripts denotes the particle index, not a tensorial index. 
The pairwise interaction is modeled by a standard 9/3 Lennard-Jones potential 


. ; 9 6 
ui el (2) -3(=) rij < Te, (2.62) 


Fij Fij 


with r;; as the distance between particle i and j. The parameters o, and eli 
describes the position and value of the minimum of the potential. Both can be 
dependent of the type of particle i and j, but for the sake of clarity indices are 
omitted. The Coulomb potential follows by 


rij < Te, (2.63) 


with the charge q; of a particle i, and the dielectric constant £]. The energy- 
conversion constant C is dependent of the particle types. Uj as well as us 
are truncated after a cutoff distance r. to reduce the computational cost, with a 
negligible error. 


Intermolecular interactions The intermolecular interaction contribution 
Uir of the potential energy function can be further decomposed by 


inter _ B A D I 
ills ae E Ut L Umit E Ue (2.64) 
ij ij, ijk, j,k, 
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(b) 


(d) 


Figure 2.6: Different contributions to the intramolecular potential energy: (a) Energy due to bonds 
stretching between atoms i and j at distance r;j, (b) due to an angle 6;;, between three 
atoms i, j and k, (c) due to a dihedral angle @; ikl of four atoms i, j, k and /, and (d) due 
to an improper angle 7; jx; of four atoms i, j, k and I. 


This involves the energy contribution due to bonds stretching between atoms 
UR, due to an angle 6; ;, between three atoms (Uf), due to a dihedral angle @; jx; 
of four atoms (Un: and due to an improper angle X;jkı of four atoms (Ula) 
An illustration of the bond length and the different angles is given in Figure 2.6. 
Note that the indices i, j,k,/ are chosen from the set of atoms that are actually 
bonded with the corresponding topology (see Figure 2.6). The individual terms 
are based on the COMPASS force field [101-103] and are completely listed in 


Appendix A.l. 
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molecular dynamics! 


In the context of thermoset polymers, early studies on epoxy systems were 
conducted by Barton et al. [104, 105], and also material property prediction of 
thermoset polymers, based on MD simulations, were performed [106-109]. To 
be able to model the reaction of thermosets, force-fields like REAXFF [110] or 
RMDff [111] were invented. Despite the ability to model reactive processes, 
the increase of the computational cost hinders the investigation of bigger sys- 
tems with more realistic number of molecules. In contrast to this, custom 
scripts were often developed in order to generate systems in a preprocessing 
step [112, 113]. Furthermore, empirical modeling of reactive processes in clas- 
sical molecular dynamics were introduced [76, 114-118]. These approaches 
compare the pre-reaction topology, and if the reaction occurs, the topology is 
updated according to a post-reaction template. Recently the REACTER frame- 
work [116-118] enabled multiple reactions during a continuously running sim- 
ulation, based on, e.g. a distance, orientation, user-specified, or a more ad- 
vanced Arrhenius type criterion. This allows massive, parallel simulation of 
thermoset polymerization. Schwab and Denniston [76] develop a similar ap- 
proach to model the polymerization of a UPPH resin system, using an Arrhe- 
nius type criterion. They were able to investigate the resin system during the 
copolymerization process with reasonable computational effort, and determine 
effective properties. 


'This section is based on the work of Schöller et al. [2]. Minor linguistic changes and additions 
have been made. 
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Urethane chain 

Radical chain 

Coupling agent 

Surface atom 

Urethane crosslink 
Connection to coupling agent 


Connection to urethane 


Fiber surface 


Figure 3.1: Schematic interface between glass fiber surface, fiber size and the resin. The legend 
on the right indicates the individual components. 


3.1 Constituents 


A schematic showing how the constituents are arranged in the system is given 
in Fig. 3.1. The surface of the glass fiber forms a passive substrate at the 
bottom for the sizing layer which is then followed by the bulk UPPH resin at 
the top. These constituents are described in more detail below. The goal is 
to construct a slightly simplified, but reasonably realistic, chemically bonded 
structural model of the full system. 


3.1.1 Glass fiber surface 


The material class for fibers used in FRPs can vary widely. One of the most 
commonly used is silica based glass, which consist in general of a complex 
chemical structure [119]. In addition, the composition of a glass surface can 
differ greatly from the bulk of the glass fiber [64, 70]. An E-glass is assumed 
for this work as it is commonly used in FRPs and is chemically simpler than 
other choices [120]. The total composition of all components on such a surface 
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Element Atomic composition (%) Atomic weight 


Si 22.5 28.086 
(0) 59.5 15.999 
Ca 5.5 40.078 
Al 6.9 26.982 
Mg 1.0 24.304 
Na 0.6 22.990 
B 4.0 10.806 


Table 3.1: Components of the glass fiber surface [70, 121]. 


has been measured by XPS scans performed by Liu et al. [70], and listed in 
Table 3.1. 


The connection between the fiber surface and the sizing agents is primarily 
established through bonding with silicon atoms. The other surface atoms play 
an otherwise passive role in this system. Hence, as a further simplification, the 
fiber surface is represented solely by its silicon atoms. Forming essentially a 
static substrate for the system, a layer of silicon atoms needs to be generated 
that provides a reasonable spacing and configuration to account for the whole 
composition of the fiber surface. 


3.1.2 Resin system 


This section introduces the resin system used in this work. The reader is re- 
ferred to the work of Schwab and Denniston [76] and Verleg [122] for a more 
extensive discussion of the UPPH resin. In order to provide the resin with 
the functionality to perform the copolymerization, it consists of three main 
components: Isocyanates, unsaturated polyester, and peroxide molecules. For 
the latter TRIGONOX® C is assumed, where only relatively few molecules are 
present, since it serves only as an initiator for the radical reaction. The iso- 
cyanate LUPRANATE® M 20 R consists of methylene diphenyl thiocyanate 
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(MDI) and their polymeric variants (P-MDI), with different functionality. Fi- 
nally, the unsaturated polyester DARON® AQR 1009 consists of molecules of 
different functionality diluted by styrene. For a complete composition with 
functionalities and molecular weights the reader is referred to Table 3.2 and 
Schwab and Denniston [76]. In the following, the two-stage reaction is de- 
scribed as used in a typical sheet molding compound (SMC) process, in which 
glass fiber-reinforced UPPH is employed [7]. 


Polyurethane reaction On a carrier film, the mixed resin is applied. In 
addition, chopped glass fibers are added. At room temperature, the hydroxyl 
groups of the UP begin to react with the isocyantes groups of (P)-MDI, creating 
crosslinks between these components (cf. Figure 3.1). This polyaddition leads 
to long polyurethane chains for which the material changes from a fluid to a 
highly viscous, more rigid state. This B-staged material can be more easily 
transported, stored and cut for the final curing. 


Radical polymerization The B-stage material is placed in a press with pre- 
heated molds. Due to the increase in temperature and pressure during compres- 
sion, the peroxide group begins to cleave, producing radical oxygen atoms. The 
free radicals attack the double-bonded carbon of the styrene and polyurethane 
molecules. By creating a more stable bond with the carbon, the other carbon of 
the double bond itself becomes a radical. Thus, a polymer chain begins to form 
and propagates through the material, forming a network of polystyrene that 
crosslinks the polyurethane chains formed by the first reaction. This results 
in an interpenetrating polymer network (IPN) of polyurethane and polyester 
chains, as shown in Figure 3.1, which gives the material its final strength and 
properties. In the system investigated in this work, the selected coupling agent 
(cf. Section 3.1.3) of the size also takes part in this reaction. This leads to 
increased adhesive bonding, as strong covalent bonds ultimately connect all 
components, including the fiber surface. The chains are terminated by the com- 
bination of radicals at the end of two chains. At the macroscopic level, the 
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Figure 3.2: Skeletal formula of the hydrolyzed y-MPS. 


material begins to cure due to heat and pressure as it flows to form the final 
component. 


3.1.3 Fiber sizing 


For the selection of glass fiber sizing there exists a large number of possible 
components and compositions. Thomason [62] describes sizing as a black- 
box technology, since size formulations are kept secret by suppliers and the 
understanding of glass fiber sizing in the literature is quite limited, fragmented, 
and usually only discussed generically. Because the size must perform a variety 
of tasks, a size used in the industry may consist of ten or more components. 
Thus, components such as a coupling agent, film former, lubricant, emulsifier, 
and other processing aids, among others could be used [123]. However, most 
of them are present only in relatively low concentrations [65]. 


The main components are therefore the film former and the coupling agent. 
The coupling agent improves the adhesion of fiber and polymer matrix. There- 
fore, it must react with the fiber surface as well as providing a functional group 
that reacts with the resin. The film former is usually a polymer that is mainly 
intended to protect the fibers during processing in the FRP production pro- 
cess [124]. The film former is usually chosen to be compatible with, if not 
identical to, the matrix material used. 


In this work, the complex sizing system is simplified using the following as- 
sumptions: 
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e The sizing components are reduced to film former and coupling agent. 


e Partly due to lack of information to the contrary, and the variety of other 
formulations possible, it is assumed that the film former is the same poly- 
mer as the matrix resin, hence UPPH, cf. Section 3.1.2. 


Only acommon coupling agent compatible with the resin is chosen, mak- 
ing our sizing somewhat generic. Most actual sizing would consist of a 
blend of different coupling agents, but these vary amongst different man- 
ufacturers and their identity and concentration are typically proprietary. 
Thus, it would be difficult to include these components in a model which 
would fit any particular system more realistically. 


3.1.3.1 Coupling agent 


Most coupling agent for glass fibers are organofunctional silanes, and their 
general structure consists of R'—Si(OR7)3, where R! provides the ability to 
react with the matrix material, creating mainly strong bonded links between 
matrix and the coupling agent, cf. Figure 3.1. R? is usually a methyl or ethyl 
group [62], and their hydrolyzed version is able to bind to the fiber surface. 
Together with the bond to the matrix, this improves the overall adhesion of the 
polymer matrix and glass fiber. 


Thomason [62] extensively studied coupling agents being used based on lit- 
erature and patents. From the manifold number of different available silane 
molecules, he concludes that the industry appears to have focused mainly on a 
few silanes where R! contains either an amino, epoxy, methacryloxy or vinyl 
functional group. The most common silane is y-aminopropyltriethoxysilane, 
which is normally used for thermoplastics and sometimes also for polyester 
and epoxy polymers. In contrast, the primary coupling agent for polyester 
appears to be y-methacryloxypropyltrimethoxysilane (y-MPS). Furthermore, 
y-glycidoxypropyltrimethoxysilane is used for epoxy and multicompatible 
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Figure 3.3: Hydrogen bonds of the coupling agent y-MPS at the glass fiber surface (left) and with 
other coupling agent molecules (right). 


polymers [62]. As in this work a UPPH resin is considered (cf. Section 3.1.2), 
the unsaturated polyester compatible coupling agent y-MPS is chosen. 


3.1.3.2 Application of the size 


In order to apply the size to the glass fibers, it needs to be dissolved, usually 
using water as the solvent. In this process, y-MPS of the size is hydrolyzed 
to dissolve the methyl groups and yield methanol and R!—Si(OH)3, with the 
complete formula shown in Figure 3.2. The aqueous size is then applied to 
the hot glass fiber by rollers, shortly after the liquid glass has been extruded 
through bushings to form the glass fiber [123]. The applied size then dries on 
the glass fibers, and the hydrolyzed y-MPS undergoes a condensation process. 
In this process, the hydroxyl groups form hydrogen bonds with other hydroxyl 
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groups. Potential reaction partners are Si-OH groups on the silicate surface of 
the fibers, but also other coupling molecules. In a subsequent step, a covalent 
bond forms while they lose a water molecule. Figure 3.3 shows the bonding of 
the y-MPS to the fiber surface and other coupling agents. Ultimately, a small 
amount of the coupling agent becomes strongly bonded to the surface, while 
others form oligomers. 


3.2 Molecular modelling description 


The models and subsequent simulations described in this section are all imple- 
mented in the open-source molecular dynamics software LAMMPS [125, 126]. 
The CoMPASS force field [101-103] is used to model the molecules and their 
interaction. In general, the potential of the force field consists of various con- 
tributions that can be grouped into two categories: The inter-molecular interac- 
tion is taken into account via individual potentials for bonds, angles, dihedrals, 
and impropers. On the other hand, a 9-6 Lennard-Jones potential models the 
van der Waals forces, and together with a Coulombic potential represent the 
intra-molecular interactions. For a more detailed discussion of the potentials, 
the reader is referred to Sun [103]. 


3.2.1 Fiber surface 


During the main simulations, the atoms on the fiber surface will be rigid. How- 
ever, to generate the surface and provide the system with the basic layer of 
silicon atoms, the use of complex potentials, including force fields, is omitted. 
Instead, a standard 12/6 Lennard-Jones potential 


E=4e (2) -(2) |. 6.1) 
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for the Si atoms is introduced, with € = 1 kcal mol7!, and o = ¥/2.0r9. Based 
on Si-O-Si bond lengths, the equilibrium spacing for the potential is set to 
ro = 3.28 A, and a cutoff length r. = 10A. Since the main objective of this 
work is to show the interaction between sizing and resin, the influence of the 
simplification of the E-glass fiber surface is assumed to be negligible. 


Domain Component No of mol. 
functionality = 2 630 

(P-)MDI functionality = 3 72 
functionality = 4 324 

Resin 1x basic structure 486 
UP 5x basic structure 36 

10x basic structure 198 

Styrene 4716 

Peroxide 72 

y-MPS hydrolyzed 2211 
functionality = 2 444 

(P-)MDI functionality = 3 51 
functionality = 4 228 

Size 1x basic structure 343 
UP 5x basic structure 25 

10x basic structure 140 

Styrene 3324 

Peroxide 51 

Surface y-MPS attached & hydrolyzed 384 
Si 4890 


Table 3.2: Composition of the different layer of the domain. Note that the number of size mole- 
cules are not the final ones, as the size layer is processing during the generation. 
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3.2.2 Reaction algorithms 


Three different reaction processes are examined in this work: First, a condensa- 
tion reaction in which the hydrolyzed coupling agent forms oligomers. Second, 
the polyurethane addition of the resin results in the B-stage material. Further, 
the final curing occurs by a radical polymerization, which crosslinks the com- 
ponents, including the coupling agent, and thus the fiber surface. Because 
different reaction algorithms are used for the condensation and polymerization 
reactions, they are introduced separately in the following. 


Condensation To generate the condensation reaction, the REACTER frame- 
work [116, 117] is used. Based on pre- and post-reaction templates, this frame- 
work offers great flexibility in modeling reactions. Although it is possible to 
incorporate various conditions, such as an Arrhenius type condition, only a 
simple distance criterion is used for the condensation: If the initiator atoms are 
within a certain distance, a drawn random number is compared with a given 
probability. If the reaction takes place, the topology is updated according to the 
specified reaction templates. Since the main goal of this reaction is to model 
the formation of coupling agent oligomers from monomers, this approach is 
considered sufficient. In addition, REACTER allows the removal of atoms from 
the system during the reaction. This allows easy treatment of the water, which 
is a product of condensation that would otherwise have to be removed by other 
cumbersome techniques (such as via a diffusion-driven drying process). 


Polymerization Following Schwab and Denniston [76] a more involved 
approach is used for both polymerization reactions. Similar to the REAC- 
TOR framework, potential reaction sites and post-reaction topology are defined 
based upon templates. Instead of only a simple distance criterion, an Arrhenius 


F = exp (5) (3.2) 


type equation 


3.3 Results 


with an activation energy AE, the Boltzmann constant kg and the temperature 
©, is considered. The activation energy AE describes the difference between 
the energy of the bond and the energy at which a bond dissociates, or alter- 
natively a bond is formed. For bond dissociation, the F must be less than a 
number drawn at random between 0 and 1. Conversely, F must be greater than 
the random number to form a bond. In addition to this basic algorithm, Schwab 
and Denniston adapted the approach to model the polymerization of UPPH ina 
more appropriate way: In the polyurethane reaction, a hydrogen atom is trans- 
ferred from a hydroxyl group to a nitrogen atom. To avoid instabilities, they 
added a transient bond that excluded these atoms from all potentials except the 
Coulomb potential, when the main bond between carbon and oxygen atoms 
is formed. This leads to a smooth transition of the hydrogen atom towards 
the nitrogen atom until the bond changes into a standard one with the full po- 
tential. Furthermore, Schwab and Denniston introduced additional Coulomb 
forces to the initiator atoms. The existence of an interaction between these 
atoms is physically motivated by these sites being electrophile or nucleophile 
centers. However, the form of the interaction (Coulomb) is artificial. These 
artificial charges accelerate the polymerization and allow it to be simulated on 
the timescales that a classical MD simulation can handle, and is therefore also 
used in this work. The charges are chosen small enough so that properties of 
the resulting bonded network does not appear to be altered by the accelerated 
dynamics [76]. 


3.3 Results 


As in the actual system, the manufacturing of the fiber, application of sizing, 
and embedding in the resin are all done in separate processes. In the following 
sections, the generation of the fiber-resin system is described in successive 
steps. The setup and conditions of each step are separately introduced, and 
their results are discussed: 


45 


3 Two-stage polymerization using molecular dynamics 


> 


SUR 


x 


temporal evolution 


> 


temporal evolution 


46 


E Pure resin | ® Size | E Fiber surface | E Cross-links | E Radical chains 


The partition of the domain in the individual layers in the left, the system after the polyurethane reaction (50% conversion rate) 
in the center, and after the radical reaction in the right. The individual components are highlighted in the corresponding color. 


In the top, sections of the domains are shown at successive time points in the reactions. 


Figure 3.4 


3.3 Results 


System Step Time step (fs) No. of steps 


Resin COMpESSSIOn 10 200000 
Equilibrating 200000 
Compression 100000 
Sizing Equilibrating 1.0 200000 
Condensation 1000000 
i 1 
Gombang re 10 00000 
Equilibrating 200000 
Polyurethane reaction 2-5 x 107? =æ 5700000 
Radical polymerization 5x107? =~ 5140000 


Table 3.3: Number of steps and time steps for each simulation step and the different systems. For 


the polyurethane reaction, the time step is increased linearly up to a conversion level of 
50%. 


Two rigid glass fiber surfaces, consisting of silicon atoms, representing 
the complex structure of a real E-glass fiber, are generated. 


The sizing layer, consisting of UPPH resin and coupling agent, is placed 
between these surfaces with some pre-attached y-MPS and the conden- 
sation reaction is conducted. 


Separately, a pure UPPH resin layer is generated and compressed to an 
initial configuration. 


The surface-sizing layer and the resin layer are combined: The two sur- 
faces bound the domain normal to these surfaces, followed by two sizing 
layers and a pure resin layer in the middle of the domain. Such a system 
is also schematically shown in the left of Figure 3.4. 


The two-stage curing process, the polyurethane and radical polymeriza- 
tion, takes place subsequently. 
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The PACKMOL [127] software is used for packing of the initial molecules. In 
combination with MOLTEMPLATE [128], this allows a flexible setup of complex 
MD simulations in LAMMPS [125, 126] using force fields. A velocity-Verlet 
integration scheme is used to solve the Newtonian equations of motion [97]. 
Temperature is adjusted by a Langevin thermostat [99] and the pressure is also 
set, as described in the next subsection. The result is that the simulations are 
performed in aNPT ensemble: Constant number of atoms, constant pressure, 
and constant temperature, although the set pressure and temperature may be 
different during the different steps described below, just as they would be in 
the actual manufacturing of a real composite system. Furthermore, simulation 
parameters, such as the time step and the number of time steps, for each simu- 
lation step are summarized in Table 3.3. In the equilibration steps, the pressure, 
temperature, and domain length normal to the fiber surface of the correspond- 
ing systems were observed to determine a sufficient equilibration time. 


3.3.1 Fiber surface 


To create the fiber surface, a random initial distribution of Si atoms is placed 
in a large periodic domain, and then compressed to E-glass density. This is 
performed at 293.15K with interactions based on the Lennard-Jones poten- 
tial (3.1). For this purpose, the final density is assumed to be 2.58gcm°, cf. 
Wallenberger and Bingham [120]. In addition, the silicon represents not only 
the mass but the entire composition of the E-glass fibers surface, see Table 3.1. 
From the final equilibrated system, two layers are extracted with a thickness of 
6.56 Ä, which is about twice the equilibrium bond length of a Si-O-Si bond. 
Thus, each layer consists of at least two Si atoms in the normal direction. The 
final edge length of the (square) fiber surface is 200Ä and also determines the 
domain sizing of the subsequent steps. In these, the silicon atoms are assumed 
to form a rigid layer based on the configuration created in this step. Therefore, 
the force on these atoms is averaged for the entire layer and an additional force 
representing the ambient pressure is applied. As a consequence, the pressure 
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can be controlled while enabling the layer to move while maintaining the initial 
configuration of the surface atoms. 


3.3.2 Sizing layer 


The implementation of the reaction algorithms (cf. Section 3.2) would allow 
depicting the complete condensation reaction of the coupling agent, including 
bonding to the fiber surface as well as to other coupling agent molecules, cf. 
Figure 3.3. However, in combination with the fully atomistic approach, realis- 
tic surface coverage is not possible due to the lack of a feasible timescale. In- 
stead, some silicon atoms of the inner surface layers have pre-attached coupling 
molecules before the actual condensation reaction takes place. This simplified 
approach seems legitimate, especially since the first chemically absorbed layer 
of the sizing appears to be fairly well understood [70, 129-131] and reasonable 
assumptions can be made about the spatial distribution of y-MPS. Based on the 
area occupied by a coupling agent molecule on the fiber surface, reported by 
Miller and Ishida [129] to be 0.59 nm?, a method for placing these pre-attached 
molecules was developed: The occupied area is converted into a mean separa- 
tion distance of the y-MPS molecules. Assuming a uniform distribution and a 
circular area zr’, this results in a minimum distance given by the radius r. For 
the first pre-attached molecule, a random atom is selected from the generated 
surface layer. Thereby, only the silicon atoms in the inner half of the layer are 
considered. Next, the Si atom with the least distance to all other pre-attached 
molecules is searched, keeping the minimum distance of /0.59/x nm. This 
atom is converted from a simple silicon atom into one with a pre-attached vari- 
ant. This search for atoms with minimal but valid distance is repeated until 
no further pre-attached molecule can be placed. This procedure results in a 
certain amount of already chemically absorbed y-MPS on the fiber surface, cf. 
Figure 3.1. 
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Figure 3.5: Distribution of the number of functional groups of the coupling agent after the simula- 
tion of the condensation reaction and the temporal evolution in the inset. 


In contrast, the main sizing layer consists of free hydrolyzed coupling agent 
monomers and UPPH resin acting as a film former, as they are considered to 
be only physically absorbed initially. The exact composition of the sizing can 
vary greatly in industrial applications. In this work, 30 wt% hydrolyzed y-MPS 
and 70wt% film former are assumed, which is within the range reported in 
the literature [62, 123, 132]. The sizing molecules are placed between the fiber 
surface layers with the previously attached coupling agent molecules. The com- 
position of the film former corresponds to the pure resin layer, and the detailed 
number of individual molecules of the sizing as well as the surface layer are 
listed in Table 3.2. Since the process attempts to mimic real-world processing 
of FRP, the condensation reaction takes place without the pure resin system, 
as the sizing is applied and dried during the manufacture of the glass fiber. 
In addition, if it is assumed that the resin is not involved in the condensation 
reaction, the computational effort can also be reduced. Finally, these mole- 
cules, the hydrolyzed y-MPS and the film former, with the exception of the 
pre-attached y-MPS, are again randomly positioned in a large domain and sub- 
sequently compressed to an approximate density of 1 gcm~?. This is followed 
by equilibration at 373.15 K with an ambient pressure of 1 atm applied via the 
fiber surfaces. 
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3.3.3 Condensation reaction 


In the condensation reaction, the hydrolyzed y-MPS monomers undergo a pro- 
cess to form oligomers. These oligomers can emerge from the physically 
absorbed coupling agent molecules of the sizing layer. In addition, the chem- 
ically absorbed sizing on the fiber surface, consisting of the pre-attached cou- 
pling agent molecules, also participates in the condensation process: Hydrogen 
bonds between the hydroxyl group and the silicon atoms form covalent bonds 
while water is being given of, see right in Figure 3.3. 


The reaction is conducted starting with the configuration of the previous step, 
at 373.15 K and | atm, to take into account the conditions of the glass fiber 
manufacturing. Furthermore, the domain is assumed to be periodic tangentially 
to the fiber surface. Normal to the surface, the rigid silicon atoms limit the 
domain. In addition, the initial reaction probability is artificially reduced and 
steadily increased during the course of the condensation process. This is done 
to ensure a numerically stable simulation, otherwise the very high reaction 
speed at the beginning would lead to instabilities. 


Figure 3.5 shows the final distribution of the number of functional groups in 
the condensed sizing. Although dimers and trimers are present, the amount of 
longer oligomers is relativ low, and most of the coupling agent is still present as 
a monomer. The inset in Figure 3.5 illustrates the development of these groups 
over time. At the beginning of the simulation, despite the lower probability, 
the reaction exhibits a high condensation rate, thereafter the rate decreases un- 
til the end of the simulation. Although it was observed that the number of 
oligomers increases with time and their functionality also increases, the distri- 
bution of the coupling agent does not completely converge to a steady state (i.e. 
the lines in the inset of Fig. 3.5 are not perfectly flat at long times). Several 
reasons seem plausible for the high proportion of monomers in the final distri- 
bution: This can be attributed to the basic reaction algorithm, e.g. the water is 
removed instantly, whereas technically the water has to diffuse though the sys- 
tem to a free surface. As the condensation is a reversible process, these water 
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molecules could undergo multiple reactions, influencing the final distribution 
of y-MPS. Moreover, the full atomic approach used in this work does not al- 
low for the very long timescale for such a condensation reaction in real life. A 
possible solution could be to omit the condensation reaction and instead work 
with oligomers from the beginning. However, since detailed information on 
the distribution of coupling agent oligomers, especially in realistic sizes, is not 
available in the literature, this approach would raise further issues. In contrast, 
a study of the condensation stage, e.g., based on a coarse-grained approach 
such as a united-atom model, could potentially go to longer time scales and 
provide detailed information about such a distribution. However, such a model 
would still need to be parameterised based on the fully atomistic model used 
here. Although such a study is highly desirable, it would be a non-trivial under- 
taking and is beyond the scope of this work. Moreover, it is expected that other 
assumptions in the introduction of the model lead to a higher uncertainty in the 
results. Therefore, the limitations of the presented approach for the conden- 
sation reaction are acknowledged, but the possibility to study the interaction 
between surface sizing and resin is nevertheless greatly extended. 


3.3.4 Pure resin layer 


In order to include a pure UPPH resin layer, a (separate) realistic block of resin 
must be created first. Therefore, the UPPH resin molecules, which are multi- 
ples of the system introduced by Schwab and Denniston [76], see Table 3.2, are 
randomly placed. As in the previous steps, the system is then compressed to an 
approximate density of 1.109gcm°, cf. [76], and subsequently equilibrated at 
293.15 K. During this process, the approach of Berendsen et al. [100] is used 
to control the ambient pressure of 1 atm, which could be used as the domain is 
assumed to be periodic in all directions. 
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Figure 3.6: Evolution of the number average molar mass M,„ and mass average molar mass Mw 
for the polyurethane reaction (a) and radical polymerization (b), over the conversion 
degree. 


3.3.5 Complete system 


In a final step, the complete system is assembled with all three layers: First, the 
surface sizing system is partitioned to form the outer layers. The final sizing 
layer is chosen to be about 50.0Ä thick. Thus, all molecules for which the 
distance of their center of geometry to the fiber surfaces is below this value 
are extracted and form the final sizing layers. Therefore, the components of 
the sizing layer differ from the initial surface-sizing system described in Table 
3.2, while approximately maintaining the relative composition. Moreover, the 
molecules of the resin layer are unwrapped in the normal direction, since the 
assumption of periodicity is no longer valid. Subsequently, the resin layer is 
placed between the both surface-sizing layers with a slight initial distance. This 
system is again compressed to an approximate density of 1gcm°, followed by 
equilibration at ambient pressure of 1 atm and 293.15 K. 


The system consists of two rigid layers of silicon atoms, followed by the sizing 
layers, with a ~% 130A thick layer of resin in the center. The resulting total 
domain has a base area of 200 x 200A and a ~ 243A large extent in the fiber 
normal direction. Based on this system, cf. left of Figure 3.4, the two stage 
polymerization reaction is carried out. 
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Figure 3.8: (a) Spatial distribution of the crosslinking density of the polyurethane reaction at dif- 
ferent conversion degrees based on a smooth representation of the data (KDE). (b) 
Spatial distribution of the relative frequency of radical chains at the end of the radical 
polymerization. 


3.3.6 Polyurethane reaction 


During the reaction, the unsaturated polyurethane reacts with the isocyante 
groups of the (P)-MDI, resulting in a polyaddition reaction that finally yields 
long polyurethane molecules. This reaction is carried out at 353.15 K and ambi- 


ent pressure [7]. In order to speed up the simulation, +1.5 charges are added to 
the initiator atoms of the reaction. As discussed by Schwab and Denniston [76], 
without auxiliary charges, the polymerization reaction would be too slow for 
the time scale used in MD simulations. Furthermore, they demonstrated the 
negligible influence of the charges on the results. 


In the center of Figure 3.4, the whole system is shown at 50 % conversion. 
In addition, the cross-links of the isocyanates and the unsaturated polyester 
are highlighted, and cutouts illustrate the development over time, i.e. the in- 
crease in cross-link density. Furthermore, the reader is referred to the ESI of [2] 
for a complete animation of the reaction. The evolution of the molar masses 
weighted by the number M, and the molar mass M,, versus the conversion de- 
gree is shown in Figure 3.6a. The detailed temporal evolution of the conversion 
degree over the simulation is plotted in Figure 3.7a and, in addition to the av- 
erage value, the median and some quantiles are also displayed. The initial fast 
reaction rate seems to converge quite quickly. Also, despite the average and 
mean values being nearly the same, it could be observed that the lowest 5 % of 
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the bins had much lower conversion degrees than the average. The domain is 
then binned based on the distance to the fiber surface to produce Figure 3.7c 
which shows the spatial distribution of the binning of the conversion degree 
over the distance to the glass fiber surface. In addition, the average value and 
a Kernel Density Estimator (KDE) are displayed. The KDE is an estimate for 
the underlying probability density function of arandom variable, and therefore 
provides a smooth function estimate of the distribution. The figure exhibits the 
highest conversion degree in the resin layer, far away from the surface. And it 
decreases steadily towards the glass fiber surface. Since the conversion degree 
refers to the relative number of reactions that took place in each bin, it is not 
possible to discuss the spatial distribution of the total number of reactions that 
occurred based on such a quantity. Therefore, Figure 3.8a shows the spatial 
distribution of cross-link density, i.e. a smooth representation of the number of 
cross-links per volume in each bin during the reaction. 


In the technical SMC process, the B-stage material represents only a semi- 
finished product, which can also be assumed to be not fully cured. In this work, 
a conversion degree of 50 % of the PU reaction is assumed as the starting point 
for the radical reaction. Especially since the definition of the conversion degree 
can vary: In this work, the number of actual reactions relative to the theoretical 
number of possible reactions is used. In contrast, an experimentally determined 
degree of conversion may differ, since the total number of possible reactions is 
generally unknown. In Figure 3.6a a drastic increase of M, at around 30-40 % 
conversion can be observed. This jump over orders of magnitude indicates the 
phase-like transition from liquid resin to rubbery B-stage material. As the ma- 
turing of B-stage is optional [7], the assumed conversion degree of 50% as a 
basis for radical polymerization is in the range of higher M,, and is considered 
a reasonable choice. A variation of this choice in combination with the inves- 
tigation of the influence on final properties is high desirable, but beyond the 
scope of this work. 


Regarding the crosslinking density, cf. Figure 3.8a, clearly different regimes 
can be identified for the sizing and resin layers during polymerization. This 
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is most likely due to the fact that the coupling agent in the sizing layer dilutes 
the UPPH resin of the sizing and thus, the final density of (P-)MDI is reduced. 
Moreover, near the fiber, the amount of y-MPS is high due to the pre-attached 
coupling agent, therefore the conversion degree is the lowest, which leads also 
to the low conversion degree in the 5 % quantile in Figure 3.7a. 


The spatial distribution of the benzene on the (P-)MDI and the unsaturated 
polyester during the polyurethane reaction is shown in the appendix A.2. While 
the distribution of the unsaturated polyester is quite smooth, the distribution of 
(P-)MDI is uneven. This behavior is less pronounced in the resulting crosslink 
density. Finally, an almost smooth transition from the resin layer to the surface 
and to the fiber of the conversion degree can be observed in Figure 3.7c. Since 
the (P-)MDI consist of the largest molecules in the system, the uneven distribu- 
tion most likely results from the generation of the system in combination with 
the limited system size. A reason why a smooth conversion degree neverthe- 
less occurs is not evident to the authors, but could be the subject of subsequent 
work. However, an experimental investigation of this very thin boundary layer 
would be very desirable as it would allow a comparison. 


3.3.7 Radical polymerization 


The B-staged material of the previous step is further processed in the radi- 
cal polymerization: Under high temperature and pressure, the fiber-reinforced 
UPPH resin flows in the molds and cures completely. Therefore, a tempera- 
ture of 418K, and a pressure of 50bar in addition to the ambient pressure is 
assumed to account for the conditions in the SMC process [7]. Moreover, the 
reaction is conducted starting with the 50 % conversion configuration of the PU 
reaction. 


In the right of Figure 3.4, the final system after the reaction is shown. The 
resulting radical chains are highlighted, and the cutouts show the evolution of 
these chains over time. As before, for a complete animation of the reaction, the 
reader is referred to appendix A.2. In addition, the evolution of the molecular 
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Figure 3.9: Spatial distribution of the radical chains length between PU crosslinks at different 
conversion degrees based on a smooth representation of the data (KDE) during the 
radical polymerization. 


weights are shown in Figure 3.6b. For the molecular weight M,, there is no 
sharp increase during the reaction, since there is no phase transition, compared 
to the PU reaction. Rather a smooth transition from rubbery to solid during the 
flow in the mold can be observed. In contrast, M„ exhibits a strongly nonlinear 
increase when styrene monomers, initially high in number but low in molecular 
weight, forms chains during polymerization. 


For the visualization of the temporal evolution of the conversion degree over 
the simulation in Figure 3.7b, the same approach as before is used. In addition, 
the vertical dotted lines mark an increase in the additional charges to accelerate 
the reaction and avoid excessive computational effort. For the start, charges of 


+1.0 were chosen, which are increased by +0.25 for each increment. Since 
the value of these charges has no direct physical effect other than hastening the 
reaction, no influence of the increase is expected if they are chosen within an 
acceptable range, cf. [76], which is not exceeded in this work. Furthermore, 
the spatial distribution of the conversion degree is shown in Figure 3.7d along 
with the average and a KDE, with the latter providing a smooth interpolation. 
Moreover, Figure 3.8b displays the relative frequency of radical chains in rela- 
tion to the distance to the fiber surface. Finally, the length of the radical chains 
and the spatial dependency, is plotted in Figure 3.9 during the polymerization. 
The length of a chain is measured by the number of carbon atoms between 
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crosslinks of the radical chains with the PU. In the resin layer, the average 
length of the radical chains is about the same. In contrast, these chains are 
longer in the sizing, as the reduced amount of PU in this layer provides less 
possibilities for crosslinking. Moreover, at low conversion degrees, the chain 
length is almost uniform and only becomes longer towards the surface at the 
end of the reaction. In contrast, the resin layer reaches its final chain length 
rather early in the simulation. The y-MPS, which may be part of the radical 
chain, allows these larger gaps between the PU to be bridged, but is less mo- 
bile compared to the small styrene molecules. Therefore, building these longer 
chains between the less present PU requires more time. 


It can be observed that some parts of the system have a much lower conver- 
sion degree, which distorts the displayed quantiles and also causes the differ- 
ence between the average and median, in Figure 3.7b. This also indicates that 
even with a longer simulation time, a complete conversion seem not feasible. 
Schwab and Denniston [76] were able to archive a higher final conversion de- 
gree and could avoid the need to raise the additional charges repeatedly. How- 
ever, in this work, the fiber surfaces restrict the movement of the molecules. 
This reduced diffusivity in the normal direction to the surface results in lower 
overall reactivity and thus lower conversion. Also, the typical styrene odor of 
the final component indicates that complete conversion is not achieved in in- 
dustrial applications. In addition, as before in the UP reaction, the definition 
of the conversion degree may differ: The conversion degree used in this work, 
based on the theoretical possible reactions, may overestimate an experimental 
determined degree of conversion. From this reduced mobility, the result is that 
the lowest conversion is near to the surface, and is highest in the resin layer 
(see Figure 3.7d). 


In Figure 3.10, the benzene of styrene is plotted during the radical polymeriza- 
tion. The benzenes are distinguished into monomeric (left) and reacted variants 
(right). While the unreacted benzene increases abruptly at the surface and at the 
interface between the size and the resin layer, the reacted styrene shows a more 
even distribution, increasing from the surface to the pure resin layer. Figures 


59 


3 Two-stage polymerization using molecular dynamics 


0.0030 1.0 
---- unreacted itd ` reacted 


0.0025 4 


0.0020 4 


0.0015 4 


Conversion degree 


0.0010 4 
5 


Density / mol cm 


0.0005 4-5 


0.0000 + 
0 20 40 60 80 100 0 20 40 60 80 100 120 


Distance from fiber surface / A Distance from fiber surface / A 


Figure 3.10: Spatial distribution of the benzene of styrene during radical polymerization. A 
distinction is made between monomeric, not yet reacted styrene (left) and reacted 
styrene in radical chains (right). 


for the distribution of all styrene benzenes, as well as for the benzenes of other 
constituents and the silicon of the y-MPS, can be found in appendix A.2. To- 
gether with the conversion degree (Figure 3.7d) and the frequency of the radical 
chains (Figure 3.8b), they indicate that molecules, especially small molecules 
such as styrene, are able to compensate for the different composition of the 
sizing and resin layers. So the conversion degree does not suddenly drop at 
the interface, rather the highly agile styrene appears to partially diffuse into the 
various layers, resulting in a smooth distribution of conversion. 


In the distribution of the radical chains, cf. Figure 3.8b, two things can be 
noted: Firstly, directly at the fiber surface there are almost no chains, but a 
high peak of chains near to the fiber can be observed. Secondly, there is an 
almost uniform distribution of the radical chains in the rest of the domain. The 
latter indicates that the presented system is able to generate a highly linked 
system via radical chains. This occurs despite the presence of coupling agent 
in the physically absorbed layer of the sizing, as the y-MPS can take part in 
the radical polymerization. While the peak close to the surface results from 
a layering effect of the pre-attached y-MPS: The functionality to react with 
the radicals have all approximately the same distance to the fiber surface. A 
more detailed spatial distribution of silicon atoms of the y-MPS during the 
polymerization reactions is provided in appendix A.2. In addition, the quality 
of the crosslinking may vary locally and is moreover quite anisotropic near to 
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the fiber surface. To illustrate these behaviors, a second animation of the racial 
reaction is provided in appendix A.2. In this, two detail views of sections at 
the surface are plotted: One for which radicals create a strong connected layer 
at the surface, with some connection normal to the surface. In the other, no 
radicals diffuse near to this chemical absorbed Y-MPS layer, resulting in only 
a few radical chains and therefore an only loosely cross-linked surface. 


3.4 Interim conclusion 


Based on molecular dynamics simulations, a UPPH resin system was extended 
by adding a fiber surface and sizing layer. From a number of possible coupling 
agents for the sizing, y-MPS was selected and a rigid silicon layer was cre- 
ated to mimic the glass fiber surface. Furthermore, the coupling agents were 
placed as monomers in the domain and a condensation reaction was modeled 
yielding dimers and higher oligomers in the sizing layer prior to the final two- 
stage polymerization. In an additional preliminary step, the three generated 
layers — resin, sizing, and fiber surface — were placed in the final system and 
equilibrated. Based on the work of Schwab and Denniston [76], the two-stage 
polymerization of the UPPH resin was conducted on this system. Subsequent 
evaluations of the quantities during the reactions were carried out. In addition, 
the system was also evaluated along the normal of the fiber surface, allowing 
a spatial analysis of the fiber-size-resin system. Due to a lack of information 
in the literature, a direct comparison of the results was not practicable. Nev- 
ertheless, the method provides an alternative approach to the characterization 
of the resin-glass interface. For further work, collaboration with experimental 
investigations would be highly desirable. Especially since the results show an 
almost smooth transition of the degree of conversion between the different lay- 
ers. This seems counterintuitive, since the different layers differ greatly in their 
composition. Moreover, anisotropic radical polymerization could be observed 
near the fiber surface, emphasizing the importance of the choice of constituents 
of the sizing, especially the coupling agent. Based on the final cured system, 


61 


3 Two-stage polymerization using molecular dynamics 


further investigation of the generated system could be performed. An evalua- 
tion of the material properties of the system, such as the crack resistance of the 
fiber-sizing resin interface, would be of great interest, but is beyond the scope 
of this work. Instead, the curing and possible crack formation during the curing 
on higher length scales, considering multiple fibers, will be considered in the 


following chapters. 
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Fiber-reinforced polymers (FRPs) consist of fibers and resin in which the fibers 
are embedded. The final microstructure can vary depending on the manufactur- 
ing process and materials used. In particular, the ratio of fiber length to diame- 
ter has a strong influence on the structure. The compression molding process, 
which is under consideration in this work, results in long, curved fibers, which 
form fiber bundles from many individual fibers, cf. Figure 2.3. 


The generation of microstructures of straight fibers up to a certain aspect ratio 
is e.g. described by Schneider [133]. This approach was recently extended to 
curved fibers by Schneider [134], allowing for higher aspect ratio and more 
complex structures. These approaches are based on the Sequential Addition 
and Migration (SAM) algorithm. Fibers are added to the system in multiple 
steps and iteratively migrated to reduce fiber overlap and fulfill other con- 
straints, such as a desired effective fiber orientation. 


In contrast, in this work an approach is presented that is based on molecu- 
lar dynamics and the Discrete Element Method (DEM). The idea of DEM, 
which is closely related to molecular dynamics, was proposed by Cundall and 
Strack [135] for the simulation of the motion of granular media and is a well- 
established method for describing motion of rigid bodies and an overview is 
given e.g. by Kafashan et al. [136]. The approach of the present work is re- 
stricted to transitional degrees of freedom of the particles and thus differs from 
the classical DEM. The latter often also considers rotational degrees of free- 
dom, and therefore also has to account for angular velocities and torques. In 
addition, the presented method uses a basic contact model of the DEM. How- 
ever, friction is neglected and potentials for bonds and angles between fiber 
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Figure 4.1: Interaction of segments in a single curved fiber. Center: A curved fiber with diameter 
D, length L consisting of N connected segments with length /. Top: Harmonic bond 
between consecutive segments with connection vector r;j and stiffness Kp as well as 
bond angle 6;;, and bending stiffness Ką. Bottom: Contribution due to a mismatch in 
orientation with a stiffness Ko and damping by a viscosity U. 


segments are introduced. The ability to prescribe a fiber orientation through a 
fiber orientation tensor has also been implemented in the model, based on the 
work of Schneider [133]. 


4.1 Fiber description 


A system with multiple curved fibers is considered. Each fiber is described by 
the diameter D, the length L and consists of N connected segments, each with 
alength/=L/N, cf. Figure 4.1. Note that the orientation of a fiber segment is 
not a degree of freedom, but is implicitly given by the position of the segments. 
For inner segments, the final orientation follows from the normalized average 
of the neighborhood connection vector, while the orientation of outer segments 
follows directly from the normalized connection vector, cf. Figure 4.2. 
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Figure 4.2: Connection and orientation vectors, in a fiber and between segments of different fibers. 
The vectors r;j,r jg describe the connection between the centers of the individual seg- 
ments, while the unit direction vectors dj, dz are averaged from them. The connection 
vector p; describes the shortest vector between the axis of segments of different fi- 
bers. 


Potential energy function Analogous to classical molecular dynamics, cf. 
Section 2.5, the potential energy function U can be introduced by 


U = LU" +E goms + £ u 53 ur, (4.1) 


i>j i>j>k i>j 


with a term U” due to overlapping segments, due to energy stored in the bond 
of sequential connected segments (URA), and due to bending of segments an- 
gle Ua), In addition, an energy associated with a mismatch in fiber orien- 
tation with respect to a prescribed orientation is added via U; = Note that the 
subscript indices ii denote the interaction between segment i of one fiber seg- 
ment and segment i of another fiber, cf. Figure 4.2, while the indices ij and ijk 
denote two or three consecutive segments of a fiber. The corresponding sums 
are performed over sets so that all possible interactions of different fibers as 
well as interactions between individual fibers are considered. 
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A basic contact model for the potential of overlapping fiber segments is given 
by 


over __ 
Ug = 


1 ;-D), Ip; <D 
iy a Ys [pal as 


2 lo, else, 


with the shortest connection vector p;; between the axis of segments of different 
fibers i and 7, the norm of a vector |-|, and a stiffness K, cf. Figure 4.1 and 
Figure 4.2. For the computation of the shortest connection vector the algorithm 
of Pournin et al. [137] is used. The potential for energy in connected fiber 
segments is given by 

Uh = 5 Ke (|r|), (43) 


with the connection vector r;; and stiffness Kp, see Figure 4.1. This formula- 
tion penalizes any deviation from / of the length of a fiber segment. For an 
angle 6; ;, between three fiber segments i, j and k, the energy follows by 


(4.4) 


1 2 
ijk = > Radix, 


with a bending stiffness X, that accounts for the energy stored in the fiber due 
to curvature. Furthermore, the energy for the orientation of a fiber segment is 
given by 
| $ 
Ur = „Ko (A— A) 2 (r;®rij), (4.5) 


with a parameter K, and the prescribed 2™ order fiber orientation tensor A. 
The actual fiber orientation of the system follows with 


= 1 
A= Ñ E dij ® dij, (4.6) 


i>j 


which consists of a sum over the number of defined orientation vectors Ñ with 
with the corresponding indices i and j. 
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Ip; 


Figure 4.3: Implemented force due to overlapping fiber segments. In contrast to the potential 
ur , a second stiffness Ke is introduced. While Ke is applied from a distance De, the 
contribution of K is added below the distance D. 


Viscous damping Based on classical molecular dynamics, the equations 
of motion can be derived from Equations (2.49) and (2.51). Together with a 
velocity Verlet integration, cf. Section 2.5.2, an NVE ensemble follows. To 
remove energy from the system, a viscous damping force with a viscosity u is 
added and the final force follows by 


=> (= Hau). (4.7) 


4.2 Results 


The velocity Verlet integration requires an evaluation of the forces. In addi- 
tion to the force derived from Equation 4.7, an additional force has been added. 
The force due to the overlapping of the segments, cf. U", was extended by a 
second term with reduced stiffness but increased distance, cf. Figure 4.3. This 
additional force should mimic the repulsive force of e.g. a classical Lennard- 
Jones potential in order to accelerate the generation of non-overlapping struc- 
tures. 
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4 Structure generator for FRPs 


7 + 1 
0 500 1000 1500 2000 
Step 


Figure 4.4: Left: Evolution of the energies due to overlapping particles Uover, the kinetic energy 
Kin and the fiber volume fraction vr during the iterations of the generation procedure. 
Right: Generated square 2D microstructure with unidirectional fibers and a size of 
200 um and 50 % fiber volume fraction. 


Description Symbol Value 
Fiber diameter ai 
De 15 um 
K i= 
Stiffness to 
K; 1x10? $ 
Viscosity u 0.5 kg 
Time step At ls 


Table 4.1: Parameters for the generation fiber structures. 


In this work, three different simplifications of a real fiber-reinforced micro- 
structure are considered and described in the following sections. Furthermore, 
all domain boundaries are considered to be periodic in all directions during the 
generation. 
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4.2 Results 


Two-dimensional unidirectional fibers Assuming a matrix reinforced 
by straight unidirectional fibers with a high aspect ratio, the fiber descrip- 
tion can be reduced to the cross sections of the fibers. The generally three- 
dimensional problem can be reduced in two dimensions. Furthermore, the fiber 
description does not need to consider multiple fiber segments, angles or orienta- 
tions. Ultimately, the problem of generating non-overlapping fibers is reduced 
to placing non-overlapping circles. Therefore, the introduced potentials and 
forces are simplified accordingly. To avoid inappropriate initial placement of 
the fibers, the fibers are randomly placed in a much larger domain. In the first 
step a compression of the domain is performed. In a subsequent step, the sys- 
tem is equilibrated in the sense of a classical MD simulation, which leads to the 
final microstructures. For a complete list of parameters, the reader is referred 
to Table 4.1. 


Figure 4.4 shows an example of a generated microstructure. The right part 
of the figure shows the final cross section of such a two-dimensional unidirec- 
tional fiber structure with a volume fraction of 50 %. The left panel depicts the 
evolution of the energies due to overlapping particles Uover, the kinetic energy 
Kxin and the fiber volume fraction vr during the iterations of the generation pro- 
cedure. The final value of Uover indicates that the fibers do not overlap and the 
remaining energy can be associated with the force due to the low stiffness Ke, 
cf. Figure 4.3, leading to a non-zero kinetic energy even for long simulation 
times. 


Two-dimensional curved fibers As an alternative to unidirectional fi- 
bers, the fibers can be assumed to be placed in one plane. This allows account- 
ing for curved fibers, along with energies in bonds, angles, as well as different 
fiber orientations. Still it allows a reduction of the general system to a 2D sys- 
tem. Therefore, the number of possible interactions is greatly reduced and the 
individual calculation can also be simplified. The parameters for generating 
the structures are listed in Table 4.1 and Table 4.2. In addition, an orientation 
is chosen with 90 % probability in the preferred direction. 
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4 Structure generator for FRPs 
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Figure 4.5: Evolution of the different energies and the fiber volume fraction vr during the iterations 
of the generation procedure. 


Description Symbol Value 
Fiber length L 200 um 
No. segments N 20 (2D), 10 (3D) 
Ky 1.0 = 
Stiffness Ka 1.0 x 10! kaum! 
Ko 1.0 x10! $ 


Table 4.2: Additional parameters for the generation of curved fiber structures. The remaining 
parameters are listed in Table 4.1. 
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4.2 Results 


Figure 4.6: Generated quadratic 2D microstructures with curved fiber. Smaller 200 um systems 
with 25 % (a), 35% (b), and 50% (c) fiber volume fraction. And a larger 600 um 
system with 25 % fiber volume fraction (d). 


Figure 4.6 depicts generated two-dimensional curved fiber structures for dif- 
ferent domain sizes and volume fractions. Figure 4.5 shows the evolution of 
the energy contributions as well as the volume fraction over the simulation 
time. The generation procedure consists of three steps. The fibers are ran- 
domly placed in a large system to avoid overlapping fibers, which can lead to 
unstable behavior. After compressing the domain, the system is equilibrated to 
obtain the final microstructure. Although the two-dimensional final structures 
represent fibers with a circular cross section, they are treated differently in this 
work. The reduction to the two-dimensional case is chosen mainly to reduce 
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4 Structure generator for FRPs 


Figure 4.7: Generated 3D microstructures with curved fiber and a system size of 360 x 120 x 
120 um and 25 % fiber volume fraction. 


the computational effort in subsequent steps. As a disadvantage, the proper- 
ties of round fibers are lost and instead an infinitely thick structure is depicted 
based on the generated structures. 


Three-dimensional curved fibers Finally, the introduced fiber descrip- 
tion can be considered in the general three-dimensional case. For the latter 
case, the parameters for generating the structures are listed in Table 4.1 and 
Table 4.2. Note that the number of segments is reduced to 10 to limit the 
number of degrees of freedom compared to the last section. Moreover, the 
orientation in the preferred direction is chosen to 80 % and 10 % in both other 
directions. In addition, the size of the simulation domain is modified accord- 
ing to the preferred direction and the aspect ratio of the fiber. An exemplary 
final three-dimensional fiber-reinforced structure with curved fibers is shown 
in Figure 4.7. 
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5 Curing simulations of FRTS on 
a micro scale 


After considering the curing of UPPH resin in combination with a glass fiber 
based on molecular dynamics simulation in Chapter 3 at the nanoscale, the cur- 
ing simulations are performed at the microscopic length scale. First, based on 
continuum mechanics, the heat equation is derived and a curing model, along 
with curing shrinkage and thermal expansion strain, is introduced. Finally, var- 
ious fiber-reinforced volume elements (see Chapter 4) are generated to study 
the curing degree, temperature, and stress distribution during a compression 
molding process. 


5.1 Model 


The curing model of this work is based on the work of Schwab [138]. However, 
instead of introducing the phasefield and the degree of cure as separate degrees 
of freedom, a derivation based on a thermomechanically weakly coupled theory 
is chosen, cf. e.g. Prahs et al. [139]. Based on such a model, a source term 
based on the curing model and a strain decomposition are introduced to account 
for thermal strains and curing shrinkage. 


73 


5 Curing simulations of FRTS on a micro scale 


Clausius-Duhem inequality In the following, following the small de- 
formation theory, a constant mass density p is assumed. Furthermore, the 
Helmholtz free energy is decomposed by 


Y= Wat Yo, Wer # Wei(9), Wo # yole), (5.1) 


into a purely elastic and a purely thermal contribution. The CDI, cf. Sec- 
tion 2.2.2, follows by 


o-é—pw—pnde >0 (5.2) 


and, taking into account equation Equation (5.1) it can be rewritten as 


Oy\ . OyW\, 4:9 


This results in the potential relations 


_ oy _ oy 


cf. eg. [85]. 
Heat conduction The localized balance of internal energy, cf. Equa- 
tion (2.15), is given by 

pé=pr+oa-é-—div(q). (5.5) 


With é = W+n}0+n6 and using the derived potential relations, the balance 
can be rewritten as 


ow . Əy., yw, Ow, oy 
Poe Poe Paoa fae ae 


-é€—div(q). (5.6) 
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5.1 Model 


This can be further simplified, leading to 


ay 


P ggg? = Pr -div(a), (5.7) 


and with the definition of the specific heat capacity at constant volume cy = 
2 
=P SE 5, the heat equation is given by 


cð = pr—div(q). (5.8) 


Thermal jump conditions Analogous to the mechanical jump conditions, 
cf. Section 2.3.3, jump conditions can also be formulated accounting for ther- 
mal measures. The effective temperature gradient g and the effective heat flux 
q are decomposed by a linear interpolation with the order parameters given by 


N N 
q=} 6%", g=} 0*9", (5.9) 


with the phase-specific fields q”, g” for each phase. In addition, the Fourier 
law 


q% = —K"1g", (5.10) 
with isotropic thermal conductivity «% is introduced for each phase. This is 
followed by the thermal jump conditions 


[q]! n% =0, [g] = 4%? n , (5.11) 


which limit the heat flux in the normal direction and a scalar jump in the 
temperature gradient a*P | Based on this assumption, an effective, generally 
anisotropic, thermal conductivity can be determined as a function of the order 
parameters, their orientation, and the individual thermal conductivities of the 
phases. 
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5 Curing simulations of FRTS on a micro scale 


Strain decomposition The effective strain is decomposed into phase- 
specific strains e*, and these are further decomposed by 


N 
e=} o%e%, e“ = e9 t+eg+ey, (5.12) 
a 


with an elastic contribution ea, 


curing shrinkage. The latter can be formulated by 


a thermal strain ed, and a strain E% due to 


[04 
e% = a®1 (0 — 0%), ef = ags, (5.13) 


with a reference temperature 0% 


af at Which the thermal strain disappears, a ther- 


mal expansion coefficient &°, the total curing shrinkage X and the curing 
degree €% for a phase & denoted by the superscript. In addition, the phase- 
specific stresses and strain energy densities are reformulated in terms of the 
elastic strain, yielding 


o% =C* le -el -e3] l (5.14) 
1 
t= z(e- eg —ef) -C° |e" 25-2]. (5.15) 


Reaction enthalpy With the specific enthalpy h, the degree of cure €% of 
a phase & can be defined by 


h” (t 
ca = EO (5.16) 
with the total reaction enthalpy h%'% < 0. Further the rate of curing degree is 


constrained by (4 % > 0 to satisfy the second law of thermodynamics. An ansatz 
for the source term of the heat equation (5.8) is introduced by 


N 
pr=—-Y¥ o h, (5.17) 
a 
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5.1 Model 


Kamal-Sourour model The model for predicting the curing process of 
UPPH is based on the Kamal-Sourour kinetic model [140, 141]. This approach 
was recently applied to an epoxy resin by Bernath et al. [142]. The degree of 
cure 6% € [0,1] describes the degree of conversion of the underlying polymer- 
ization from uncured state &% = 0 to the fully cured state &* = 1 of one phase 
or subregion œ. It can also be phenomenologically related to the amount of 
enthalpy provided by the polymerization, cf. Equation (5.16). The rate of cure 
follows by the ordinary differential equation 


COGEO) = (K1 (0) +k" (0) (E9) (I-64), 618 


with the exponents m,n and the kinetic functions k®%!, k%?. The latter can be 
described by the Arrhenius equations 


ES! 


RO 


Et? 
k™1(9) = A®! exp (- ) ,  k*?(0) =A% exp (F) (5.19) 
with the prefactors A®! A%?, the activation energies E%!, E% and the univer- 
sal gas constant R. 


Glass transition temperature The glass transition temperature 0%S de- 
scribes a change of state of polymers. Below 0%® the material is often char- 
acterized by brittle material behavior, while above the glass transition tem- 
perature it often exhibits viscous behavior, which can be associated with a 
rubbery state. In the context of the curing process, there is no expected 
evolution of the degree of cure below the transition temperature. Since the 
Kamal-Sourour model does not account for such a transition, the approach 
of DiBenedetto [143] is used to describe the glass transition temperature as a 
function of the degree of cure, yielding 


Q%8(E%) — 9980 AREA 


0&8” — Q@,g0 a 1— (1 — 1%) ça ’ (5.20) 
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5 Curing simulations of FRTS on a micro scale 


Description Symbol Value 
Mass density p 2575 s 
Specific heat capacity Cy 802.4 RK 
Thermal conductivity K 1.275 x 
Bulk modulus K 46.5 GPa 
Shear modulus G 33 GPa 
Thermal expansion coefficient a 1x 10 t 


Table 5.1: Material parameters for the glass fiber. Values are taken from Schwab [138] (Table 
4.12), besides the thermal expansion coefficient. 


with the glass transition temperature 0%" for no cure and fully cured 0%8°, 
respectively, and the parameter 4%, which can be determined e.g. experimen- 
tally [142] or by molecular dynamics simulations [76]. 


5.2 Results 


Although the model was introduced with a degree of cure for all phases, i.e. 
materials, in the following only the resin undergoes a curing process and solely 
exhibits curing shrinkage. Therefore, in the following, the superscript is omit- 
ted for quantities related to cure, but always denotes the quantity of the resin 
material. 


The material properties used for the glass fiber are listed in Table 5.1 and are 
assumed to be constant over the considered temperatures. For both resin and 
glass fiber, the reference temperature 0%, is assumed to be 293 K. In contrast, 
the references to the material properties of the UPPH resin used, including cur- 
ing parameters, are listed in Table 5.2. Many of the parameters are temperature 
and/or cure dependent. This dependence is modeled by linear or bilinear inter- 
polation between the tabular data provided in the corresponding references. 
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5.2 Results 
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5 Curing simulations of FRTS on a micro scale 


The maximum absolute principal stress is defined by 
c= o' 5.21 
max, |]; (5.21) 


which is used to visualize the stress states of the microstructures during curing. 
In addition, the volume average is defined by 


1 
=, [wav (5.22) 


for a arbitrary field w. 


Boundary and inital conditions The initial temperature of the simula- 
tion domain is defined as room temperature with 293 K. While the initial cur- 
ing degree is set to 1 x 107° to avoid numerical difficulties due to the singular 
point of the model at zero curing degree. The boundaries are all considered 
periodic, except in one direction, which is associated with two sides of the 
molds, denoted by d Qm. Therefore, the normal displacement are restricted in 
this direction with 

u:-n=0, ondQm, (5.23) 


with the outer normal n of the boundary. For the heat equation, a Robin-type 
boundary condition is chosen on these boundaries. The heat flux is given by 


q:n=h(0—6.), ondQm, (5.24) 


with a heat transfer coefficient h and the environmental temperature @... In the 
following, a variation of the domain size is investigated. The total thickness 
of the component is assumed to be 2mm and the heat transfer coefficient is 
chosen as h = C/Al, with the length difference Al between the simulation 
domain and the total thickness of the component, and the coefficient C = 
1x 10°”Wm-!K“!. This approach takes into account the thermal damping 
of the component not modeled and avoids overly high heat fluxes for smaller 
simulation domains. 
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5.2 Results 


0 50 100 150 200 250 300 350 


Figure 5.1: Temporal evolution of effective curing degree (C) (solid line) and effective tempera- 
ture (0) (dashed line) during the cure process for pure resin and different system sizes 
in um. 


To model the compound molding process, the simulation is divided into two 
steps: 


e Closing the molds and heating the SMC prepeg to initialize the curing 
process. This is accomplished by setting 0» to 373 K. 


e The molds are opened and the final component is cooled to room temper- 
ature. Therefore, 0, is reset to 273 K after 3 minutes. 


5.2.1 Pure resin 


In the first step, a homogeneous domain with pure UPPH resin is considered 
to see the influence of the reinforcement by the glass fibers in the subsequently 
sections. Figure 5.1 shows the effective curing degree (¢} and the effective 
temperature (0) during the curing simulation. After reaching the initial glass 
transition temperature of 323 K, the curing starts and due to the exothermic 
reaction, the temperature increases and reaches temperatures above 6... The 
curing degree reaches the final value quite early, before the cured resin slowly 
cools down to room temperature. The different domain sizes have different 
thermal inertia and heat transfer coefficients, resulting in different heating and 
cooling kinetics due to the boundary conditions. In contrast, the evolution of 
the curing degree is only slightly affected, except for a time shift. 
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5 Curing simulations of FRTS on a micro scale 
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time / s 
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time/s 


Figure 5.2: Top: Temporal evolution of effective curing degree (&) (solid line) and effective tem- 
perature (0) (dashed line). Bottom: Temporal evolution of effective maximum princi- 
pal stress (0*) (solid line) and local maximum principal stress (dashed line). These are 
plotted during the curing process for unidirectional fiber-reinforced domains with a 
fiber volume fraction of 25 % and different system sizes in um. 


5.2.2 Unidirectional fiber-reinforced 


In the next step, simulation studies with unidirectionally fiber-reinforced do- 
mains, cf. Section 4.2, are investigated. Figure 5.2 displays the effective cur- 
ing degree, the effective temperature and the maximum principal stress (0*) 
during the curing simulation of a domain with a volume fraction of 25%. In 
contrast to the pure resin, the high conductivity and relatively low heat capacity 
of the glass fiber increase the rate of the effective temperature. This leads to an 
earlier start of the cure process and the final room temperature is also reached 
earlier. The maximum principal stress is caused by the interaction of strains 
due to thermal expansion, curing shrinkage and the different stiffness of both 
materials, which changes during the curing process. The highest effective and 
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5.2 Results 


RSs 
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Figure 5.3: Distribution of maximum principal stress o* for a square domain of 400 um length 
with a fiber volume fraction of 25 % for a unidirectional fiber-reinforced system (left) 
and a two-dimensional curved fiber-reinforced system (right) after curing simulation. 


local values were observed at room temperature after curing. This indicates 
that the final curing shrinkage is the main cause of these stresses. Figure 5.3 
shows an example of the final distribution of the maximum principal stress for 
a 25 % volume fraction domain. In the resin-rich areas, the stress is mostly ho- 
mogeneous, while the fiber cross section introduces local stress concentrations. 
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5 Curing simulations of FRTS on a micro scale 
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Figure 5.4: Top: Temporal evolution of effective curing degree (&) (solid line) and effective tem- 
perature (0) (dashed line). Bottom: Temporal evolution of effective maximum princi- 
pal stress (o*) (solid line) and local maximum principal stress (dashed line). These 
are plotted during the curing process for two-dimensional curved fiber-reinforced 
domains with 25 % fiber volume fraction and different system sizes in um. 


5.2.3 Two-dimensional curved fibers 


Figure 5.3 illustrates the final distribution of the maximum principal stress for 
a 25 % volume fraction domain with a curved fiber in two dimensions, cf. Sec- 
tion 4.2. In contrast to the unidirectional fiber-reinforced system, the distri- 
bution is more homogeneous because the curved fiber domains contain larger 
resin-rich areas for a given fiber volume fraction. However, the sharp corner 
of the curved fibers introduces locally higher stresses. Figure 5.4 shows the 
temperature, curing degree, and maximum principal stress over time. While 
the temperature and curing degree are similar to the unidirectionally fiber- 
reinforced variant, higher principal stress can be observed. In particular, some 
ensembles show a peak of the local maximum stress at the beginning of the 
cure. 
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5.2 Results 
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Figure 5.5: Top: Temporal evolution of effective curing degree (&) (solid line) and effective tem- 
perature (9) (dashed line). Bottom: Temporal evolution of effective maximum prin- 
cipal stress (0*) (solid line) and local maximum principal stress (dashed line). The 
is plotted during the curing process for three-dimensional curved fiber-reinforced 
domains and different and different system sizes in um. 


5.2.4 Three-dimensional curved fibers 


Finally curing simulations of three-dimensional curved fiber, cf. Section 4.2, 
are conducted. Figure 5.5 shows the temperature, curing degree, and maximum 
principal stress over time. While the curing and temperature evolution is simi- 
lar to the two-dimensional cases the simulation demonstrate higher maximum 
principal stresses, since the full three-dimensional feature of the curved fibers 
can be depicted correctly. The maximum principal stresses are visualized in 
Figure 5.6, whereas regions with lower values are omitted and only the regions 
of high stresses are highlighted in color. 
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5 Curing simulations of FRTS on a micro scale 


Figure 5.6: Distribution of maximum principal stress o* of a three-dimensional curved fiber- 
reinforced system with 25 % fiber volume fraction after curing simulation. Red high- 
lights maximum values, while lower values have reduced opacity to ensure clarity. 


5.2.5 Comparison 


During curing, the material parameters of UPPH materials change significantly. 
However, when different volume fractions are considered, the change in e.g. 
effective conductivity and thermal inertia is dominated by the glass fiber. The 
difference in the evolution of cure degree and temperature for different abstrac- 
tion levels of the volume element or selected ensemble seems to be small. In 
contrast, the mechanical behavior seems to be different for the three levels of 
abstraction. Therefore, a study with varying fiber volume fraction is performed 
for all three variants of the fiber-reinforced system. For each type and volume 
fraction, simulations were performed on 5 different ensembles with a size of 
400 x 400 um (2D) and 240 x 80 x 80pm (3D), respectively. Figure 5.7 shows 
the maximum principal stress over the fiber volume fraction for these cases. 
Therefore, both the effective value and the maximum value are plotted and con- 
fidence intervals are provided. Almost no scatter for the effective maximum 
principal stress can be observed despite the use of different ensembles. Only 
a larger difference between the two-dimensional and three-dimensional cases 
is present. This indicates that the different geometry of the two-dimensional 
cases does not cause any difference in the effective behavior, only the change 
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Figure 5.7: Maximum principal stress o* over fiber volume fraction for unidirectional, 2D curved 
and 3D curved fibers. The effective value as well as the maximum value are plotted. 
Based on different ensembles, a confidence interval is also provided. 


to a three-dimensional strain state leads to higher stresses. In addition, the 
value decreases for higher volume fractions because the a higher fraction of 
the domain do not exhibit curing shrinkage. 


The local maximum principal stress increases with the complexity of the mi- 
crostructure. The unidirectional fiber-reinforced systems have the lowest value 
and the three-dimensional curved fiber case the highest value. The ability to 
generate complex three-dimensional stress states increases with the complex- 
ity of the represented fiber shape. Also, with higher fiber volume fractions, the 
likelihood of unfavorable fiber positions, and thus areas of stress concentration 
increases. In particular, a strong increase in stress of the three-dimensional 
system from 20 % to 40 % volume fraction can be observed. 


It can be concluded that while simplifying the real microstructure of FRTS to 
two-dimensional structures does introduce an error in predicting stresses, it nev- 
ertheless allows investigating the basic qualitative mechanism of the mechan- 
ical behavior of glass fiber-reinforced UPPH resin at reduced computational 
costs. 
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6 Phase-field modeling of crack 
propagation in heterogeneous 
materials 


In this chapter, phase-field models for crack propagation are introduced and 
their application especially to heterogeneous systems such as FRPs is dis- 
cussed. First, a crack propagation model for a homogeneous system is pre- 
sented, which is based on simple but well-established models from the litera- 
ture. Subsequently, two models for heterogeneous systems follow by general- 
ization of the homogeneous model. The single-crack order parameter (SCOP) 
model describes fracture by a single parameter for all phases, i.e. materials. 
In contrast, multiple parameters are introduced for the multi-crack order pa- 
rameters (MCOP). Along with some numerical remarks and an extension to 
a single-obstacle potential, numerical results are discussed. A qualitative and 
quantitative comparison of the SCOP and MCOP is performed. Furthermore, 
the basic Voigt-Taylor scheme is compared to the mechanical jump condition 
framework, cf. Section 2.3 and Schneider et al. [60]. Finally, the application 
of these models to FRPs will be demonstrated, restricted to systems with ex- 
ternal loads, excluding the curing process. This chapter is mainly based on the 
work of Schöller et al. [1, 3], denoted by footnotes, but also includes several 


extensions of these works. 
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6 Phase-field modeling of crack propagation in heterogeneous materials 


(a) (b) 


Figure 6.1: Schematic homogeneous domain with a sharp crack interface (a) and a diffuse crack 
interface in yellow (b). 


6.1 Phase-field crack propagation models! 


6.1.1 Classical homogeneous model 


For a material body Q € R” ne {1,2,3}, in an Euclidean space, the displace- 
ment vector u relates the position vector of a material point in the current 
configuration relative to the reference configuration. The boundary of the body 
IQ € R”! consists of a subset o Qp, on which a Dirichlet boundary condition 
is applied, by prescribing the displacement vector u = %, while a Neumann 
boundary condition is imposed on 0Qn, for which the stress vector t = t is 
specified. For both boundaries, OQp Ud Qn = OQ and dQDNAQN = O must 
apply. In addition, the body contains a sharp crack surface Se. 


'The content of this section has been taken directly from Schöller et al. [1] with minor linguistic 
changes. 
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6.1 Phase-field crack propagation models 


Free energy functional For a homogeneous domain containing a sharp 
crack, cf. Figure 6.1a, the free energy functional is proposed by 


Fu] = [pfatwar+ | Geda, (6.1) 


with a strain energy density f.)(u) and a critical energy release rate Ge, which 
describes the surface energy density y for the sharp crack surface by G. = 
2y. The introduction of the crack order parameter d. € [0,1] enables a smooth 
transition of the material state. For @, = 1, the material is fully broken, for 
6. = 0, it remains undamaged. Following Kuhn and Müller [28], this allows 
the free energy (6.1) to be rewritten by a volume integral 


1 1 
Flu, bc, V Ge] = I h(de) falu) + zC (cs. IV.) au a) dy, (6.2) 
a L 
Fu,de Vo) 


cf. Figure 6.1b. The strain energy density is degraded by the function h(@,) = 
(1- $e)’. For a general discussion of degradation functions, the reader is re- 
ferred to Kuhn et al. [144]. The energy of the crack is parameterized by @¢7 
and a gradient term with the spatial gradient V6. = grad(@,) and its norm 
|Voc| = /VGc-VGc- In addition, the length parameter £, defines the width 
of the diffuse interface. 


Strain energy density The elastic term of the free energy functional (6.2) 
is modeled by an elastic potential, assuming small deformations. Thus, the 
strain energy density reads 


1 
falu) = 57 (6.3) 
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6 Phase-field modeling of crack propagation in heterogeneous materials 


with the strain tensor € = sym(grad(w)), where sym(-) denotes the symmetric 
part of a second order tensor. In addition, Hooke’s law 


o =Cle], (6.4) 


is assumed, with the symmetric Cauchy stress tensor ø, and the 4" order stiff- 
ness tensor C. Moreover, body forces are neglected and a quasi-static behavior 
is subsequently considered. For an isotropic case, this constitutive equation 
simplifies to øo = Atr(e)1+2ue, with the Lamé parameters A and u, the 
second-order identity tensor 1, and the trace operator tr(-). The Lamé parame- 
ters can also be expressed by Young’s modulus E and Poisson’s ratio v, using 


Ev E 


een) ETTEV 


(6.5) 
As cracks typically do not evolve under compression loads, several approaches 
exist to split the strain energy density into a tension and compression part. By 
only degrading the tension part, this results in a more physical crack propaga- 
tion behavior. Early tension-compression splits were introduced by Amor et al. 
[145], Miehe et al. [146] and Henry and Levine [26]. More recent variants were 
proposed by Ambati et. al [34], Strobl and Seelig [147] and Storm et al. [148]. 
Despite their significance, none of these variants is considered here, in order 
to reduce the complexity of the proposed model. However, an extension to 
most of the established tension-compression splits could be applied straightfor- 
wardly. Applying Equation (6.3) and (6.4) to the functional (6.2), yields 


1 1 1 
Flugo VO) = | ANCL e+ 50 (Va +02) an 66 


Balance of linear momentum Following the approach of Kuhn et al. [28], 
the minimization of the free energy functional (6.6), with respect to the dis- 
placement yields the balance of linear momentum 
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div(h(d.)e) = o. (6.7) 


Evolution equation The evolution equation of the order parameter s can 
be described by an Allen-Cahn equation [89] 


SF Af (Ube, VG) ~ ( AF (U, de, V Ge) 
a u( Jð. aiv( IVO. )) , (6.8) 


EG, de =-—M 


with a mobility parameter M > 0, Vdc: =0 on AQ and the outer normal 
vector n on the boundary. 


Analytical crack profile For a one-dimensional stationary crack without 
mechanical loads, and thus with vanishing mechanical driving force, the evolu- 
tion equation according to (6.8) is able to reproduce the correct surface energies 
of a sharp interface and leads to the analytical profile [29] 


bc (x) = exp (- 2) ; (6.9) 


which is displayed in Figure 6.2a. 


0 i j 0 i j 
= 3&9, —Ep 0 Ep 3E be — Eh, 0 Ep, 
X Xx 


Figure 6.2: Analytical order parameter profile for crack (a) and solid (b), for a corresponding sharp 
crack and a solid-solid interface at x = 0, respectively. 
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Irreversibility condition Various approaches for the realization of an ir- 
reversibility condition for the evolution equation are listed in the literature: 
Bourdin et al. [27] used a Dirichlet boundary condition. In contrast, Miehe et 
al. [29] introduced a strain history function to realize the irreversibility of 
crack order parameter. More recent approaches use the augmented Lagrangian 
method [149], a primal-dual-active set [150], acomplementarity system [151], 
or the interior point method [152]. The model introduced above reproduces the 
correct sharp interface energy in the case of the analytical profile (6.9). Since 
mechanical loads are considered here, the strain energy density contributes to 
the evolution equation of the order parameter (6.8). This leads to a distorted 
profile, if crack healing is prevented completely, and therefore to an error in the 
represented energy, if the applied load is removed completely. As the objective 
of this work is to investigate the possibility of different models to reproduce the 
sharp interface energies, the irreversibility is realized by a Dirichlet boundary 
condition, where the evolution equation is not restricted in general, but all fully 
damaged points remain damaged, using a boundary condition @, = 1, where a 
point is assumed to be fully damaged by @, > 0.99 using a numerical threshold. 
This improves the numerical behavior and has a negligible effect on the results. 
Nevertheless, the models presented in the present work could also be realized 
with other approaches to the irreversibility condition, which are listed above. 


6.1.2 Heterogeneous single-crack order parameter 
(SCOP) phase-field model 


Solid order parameters To distinguish different regions of a heteroge- 


neous body, the order parameters @% € [0,1], Væ =1,...,N are introduced for 
N different regions, which can be arranged in the tuple 


b= {9',97,...,0%}. (6.10) 
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Figure 6.3: Schematic heterogeneous domain for the SCOP and MCOP models: The initial do- 
main is schematically drawn on the top left. The subdomains are visualized separately 
in the center, with a sharp crack interface. In the lower left, the SCOP model is shown 
with its diffuse crack interfaces in yellow, whereas in the lower right, the MCOP model 
is displayed with several diffuse crack interfaces in the subdomains 2% and QP in red 
and yellow. The diffuse solid interface ¢% (a) 8 (a) is indicated in green. 


A bulk region of phase @ is represented by @% = 1, while a diffuse interface is 
represented by 0 < @% < 1. As the order parameters can be interpreted as the 
volume fraction of the corresponding regions, the condition 


N 
ee =1 (6.11) 
a=1 


has to be fulfilled. By extending the free energy functional (6.1) with interfa- 
cial terms between the regions, a multiphase-field model can be formulated as 
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proposed by Nestler et al. [46]. Using a double obstacle phase-field potential 
results in an analytical one-dimensional profile for a binary interface: 


6%) = 3 (14sin(Fx)), “3 S05) (6.12) 


as shown by Steinbach et al. [153], which is presented in Figure 6.2b. It is char- 
acterized by the interface width /, or, analogous to the crack profile, by a length 
parameter £p, = 4/ 771,. Since no changes of @% are examined in this work, an 
extension of the free energy and a derivation of a classical multiphase-field is 
omitted. Instead, the analytical profile (6.12) is utilized to parameterize the 
domain, resulting in diffuse and smooth volumetric interfaces, cf. Figure 6.3. 


Free energy functional To account for heterogeneous materials, the do- 
main Q can be decomposed into subdomains Q% with constant material prop- 
erties, cf. Figure 6.3. Thus, the energy functional for sharp interfaces follows 
by 


N N 
F [u] =} f aww, Go da. (6.13) 


In the context of the SCOP model, and thus, in contrast to the functional (6.1), 


[01 
el 


phase-specific, denoted by the phase index œ. The energy of each subdomain 


the strain energy densities f% as well as the critical energy release rates G% are 
is given by an integral of energy densities over 2%, which can be parameter- 
ized using an indicator function and expanded to Q. Subsequently, the order 
parameters of the solid phase field can be used as a smooth approximation to 
the indicator function [154, 155]. This leads to a linear interpolation of the 
energy densities with the solid order parameters &, cf. for example Nestler et 
al. [46], yielding 


N N 
Ful = | Losi | Zoda. (6.14) 
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As for the classic homogeneous model, a crack order parameter @, is intro- 
duced to describe the damage of the domain smoothly, and the free energy of a 
heterogeneous body is given by 


N 
F (tu, be, Vor] = f h(oe)  o%F2 (6.15) 


eM= 


1 ho 
+5 o°Ge (ea val +92) dv, 


1 
fal = get et, a” =C% [e°]. (6.16) 


In general, the phase-inherent stresses and strains are unknown. To determine 
an overall material behavior, further assumptions have to be made. This is 
widely investigated in the context of phase-field modeling [56-59, 61]. Re- 
cently, Prajapati et al. [47] introduced a model that applies the homogenization 
scheme proposed by Schneider et al. [61] in the context of a phase-field fracture 
model. Nevertheless, for simplicity, a Voigt-Taylor homogenization scheme is 
used in this work, assuming equal strains in each phase 


e“ =eVa=1,...,N, (6.17) 


acknowledging the limited capabilities of this scheme [22, 58]. By apply- 
ing (6.17) and (6.16) to the functional (6.15) this yields 


1 N 
F |U, bc, VO] = f ZHODE [e])-e (6.18) 


1x 1 
+5 Lorat (ea lvoe? + 0) av 
a Epe 
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Note that the procedure to obtain the evolution equations of the crack order 
parameter and the balance of linear momentum is the same as for the classical 
phase-field crack propagation model in Section 6.1.1, but now including inter- 
polations of the phase-specific stiffnesses and critical energy release rates with 
the order parameters @%. 


6.1.3 Heterogeneous multi-crack order parameter 
(MCOP) phase-field model 


Order parameters As in the previous section, a tuple of order parameters, 


N 
p={0',0°,...,0%}, oS, (6.19) 
a=1 


is introduced to parameterize different regions of a heterogeneous body. In 
addition, a tuple consisting of separate crack order parameters 0% € [0,1], for 
each solid phase, 


he = {010-00}, (6.20) 


is introduced. Each crack order parameter ®° tracks the damage of the cor- 
responding solid region œ. For 6% = 1, for instance, the complete volume 
fraction @% is damaged, while other solid regions are not affected by #2, cf. 
Figure 6.3. 


Free energy functional The MCOP model combines the functional given 
by Equation (6.14), used as the basis for the SCOP model, with the individual 
degradation of the strain energy densities f% by means of multiple crack order 
parameters according to Equation (6.20). In addition, the critical energy release 
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rates G are also parameterized by the individual crack order parameters. This 
results in the free energy functional 


N 
F |, be,V be] = 1 Eons (6.21) 
ad sLonog (ea vaP + > L e)a» 


with Vd = {vo} ‚V®2,...,voN }. The free energy functional for each subdo- 
main Q“ can be identified by 


Flu, 2,992] = i 9° f% u, 62,092) dv (6.22) 
= fa tOO (eq WOE + 2-00) ) ar 
gee ee es 


f% (uw, Ge. Vč) 


With Equation (6.16) and, as in the previous model, the assumption of the Voigt- 
Taylor scheme (6.17), the functional for the whole domain and each subdomain 
is obtained by 


7° 1u,08,002]= f AOC ee (6.23) 


+5 568 (eq var? +, EET — (09°) av v, 


Mz 


o% ho) (C Tel) e (6.24) 


NI = 


R 
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Flu, be, Vo.) =| 
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Pie: (ea vog + + 0e) dy 
Eo 
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As in the previous models, minimizing the latter functional with respect to the 
displacements yields the balance of linear momentum and is therefore omitted 
here. 


Evolution equation For each crack order parameter an Allen-Cahn equa- 
tion, 


Ol on OF” 
Epe =—M 500 (6.25) 
ey, Ba a) 
99% OV OS 


with a mobility M“ > 0 is postulated. These evolution equations are defined on 
the subdomains Q%, whereas on the boundary d. Q% a homogeneous Neumann 
boundary condition with Vo“ -n = 0 applies. In contrast, the boundary ¥%, 
resulting from the smooth transition of Z (cf. Figure 6.3), has to be treated 
separately. Neither a classical Neumann nor a Dirichlet boundary condition is 
a reasonable choice. A Neumann boundary condition would enforce a certain 
flux across the boundary. For example, a zero flux would result in a crack propa- 
gation direction perpendicular to the boundary. A Dirichlet boundary condition, 
on the other hand, would constrain the order parameter ¢& and thus the crack 
propagation. More complex boundary conditions, such as absorbing boundary 
conditions, or Robin type boundary conditions, could constitute more physical 
boundary conditions. Since these are neither widely used nor trivial to imple- 
ment, an alternative approach is proposed: The evolution equation is extended 
to the whole domain Q, but the elastic driving force is restricted to Q%, and 
considered to vanish anywhere else. This results in a continuous calculation 
of the crack order parameters, not restricted to the inner boundary, each with 
similar terms, as in the homogeneous model, but with phase-specific quantities. 
Outside the subdomain 2%, the phasefield is continued in the sense of the expo- 
nential profile (6.9), reproducing the correct sharp interface energy in a diffuse 
context. The coupling of these different equations takes place solely through 
the interpolation of the strain energy density. In addition, an irreversibility 
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condition for each evolution equation is used: As in Section 6.1.1, each crack 
order parameter is kept damaged by means of the additional constraint 6% = 1, 
if 6° > 0.99. 


6.1.4 Comparison 


Classic homogeneous model Section 6.1.1 introduced a crack propaga- 
tion model based on established models from the literature. Many extensions 
and modifications to such a model can be found in the literature, which can im- 
prove the model for many applications [48, 50, 148]. Nevertheless, such exten- 
sions are avoided intentionally, reducing the complexity, as the model is used 
as the basis for the introduced heterogeneous SCOP and MCOP models. In 
the scope of this work, no implementation of the classical model is conducted. 
While the limitations of the classical model are acknowledged, many of the es- 
tablished modifications might also be applied to the novel MCOP model, while 
retaining the advantages of this approach. Nevertheless, the classical model 
could be extended to heterogeneous systems, e.g. based on element-wise con- 
stant material properties in the context of the FEM. On the other hand, the 
SCOP model with the limit of a sharp interface would yield the same results. 


Single-crack order parameter (SCOP) models The SCOP model, in- 
troduced in Section 6.1.2, is similar to models from the literature [22, 47, 48]. 
However, there are some differences: 


e Schneider et al. [22] introduced the order parameters 6° and the crack 
order parameter @, in the same tuple &. Therefore, @, also contributes to 
the summation condition (6.11). This results in a model where the phases 
show transitions to a common crack phase. In addition, the interpolation 
function h(6.) must then normalize with respect to I ¢%, while the 
crack energy is also considered in the evaluation equation of the phases, 
and vice versa. 


101 


6 Phase-field modeling of crack propagation in heterogeneous materials 


e Prajapati et al. [47] and Späth et al. [48] used a similar approach as 
Schneider et al. [22], but the crack evolution equation is assumed to be 
independent of the order parameters 6°, and is solely used to determine 
effective material properties. 


e In this work, the order parameters @% and the crack order parameter Qe 
are introduced separately, and the damage of all phases is represented 
by d., but without any phase transitions from these to the crack phase. 
This results in a model similar to Prajapati et al. [47] or Späth et al. [48], 
but without requiring a normalization of the interpolation function. 


e Due to the similarities to established models, the SCOP model is used as 
a reference model in the following 2D examples. 


N 
0=div (ed ote") A 


Eb 


; 1 ana 1 ` a a 
ey = ọ Ge (ar-a) +26 vo “Vdc — 


de 


N 
0=div (Encore) ; 


ya 1 1 dh(oe 
de e 2c 


Table 6.1: Comparison of the terms in the balance of linear momentum and the evolution equation 
of the SCOP and MCOP model?. A- denotes the Laplacian operator. 


2 An error in the derivation from Schöller et al. [1] has been corrected in this work. 
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SCOP vs MCOP model Both models introduce a tuple of order parameters 
for the parametrization of the domain, thereby accounting for the heterogene- 
ity of the body. The differences between both models are summarized in the 
following. 


e SCOP model: Only a single crack order parameter is considered. Thus, 
all regions are equally damaged. Both, the balance of linear momen- 
tum as well as the evolution equation of the crack order parameter are 
obtained by minimization of the functional .7 with respect to the total 
domain Q. 


e MCOP model: A tuple of crack order parameters is introduced. Thus, the 
damage of a region is described by its own order parameter, which allows 
a more advanced degradation of the strain energy. Moreover, function- 
als Z“ are introduced on subdomains 2%. As for the SCOP model, the 
balance of linear momentum is obtained by the minimization of the total 
functional with respect to the total domain. However, the evolution equa- 
tions of the crack order parameter are obtained by the minimization of 
the functionals F% with respect to the domain Q% of the corresponding 
crack order parameter. Each evolution equation also recovers the classic 
model, while maintaining a constant crack surface energy and many of 
the advantages of this model. 


Furthermore, the differences in the evolution equations and linear momentum 
balances are summarized in Table 6.1. 
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6.2 Numerical treatment? 


Influence of the crack length parameter £ Regarding phase-field 
crack models for simulations, it has been repeatedly shown that the length pa- 
rameter used to determining the width of the diffuse transition between dam- 
aged material points and undamaged material points has an effect on the sim- 
ulation results, see [28, 29, 48], for example. This is especially true when 
considering crack initiation processes, for which often &,, is treated as mate- 
rial property, cf., e.g. Tanné et al. [156]. Recently, Kumar et al. [157] dis- 
cussed crack nucleation in the context of phase-field modeling and promoted 
an alternative approach for nucleation treatment. If initial cracks exist and the 
transition width is compatible with the discretization grid and the domain size, 
the influence of & , is not significant [156, 158-160]. This difficulty was exten- 
sively studied in [161], and it was concluded that the length parameter should 
be considered as a material property that depends on the tensile strength of a 
material. Tanné et al. [156] derived a possible solution for the correct determi- 
nation of the length parameter. However, if the length parameter is considered 
as a material property, especially on small length scales, this often leads to dif- 
ficulties, as the compatibility between the length parameter, the discretization 
grid, and the domain size is no longer guaranteed. To eliminate this sensitivity, 
Wu et al. [162] introduced an approach for brittle materials. Since the presented 
model is a completely new interpretation of the regularized crack problem, the 
approach by Wu et al.’s is not considered, to reduce complexity. For a clear 
presentation of the new model, the disadvantages and the problem of the de- 
pendence of the simulation results on the length parameter are acknowledged 
but not taken into account. However, improving the model towards parame- 
ter insensitivity is nevertheless straight forward with the approach published in 
Wu et al. [162]. 


3The content of this section has been taken directly from Schiller et al. [1] with minor linguistic 
changes. 
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Algorithm 2: Staggered scheme for SCOP Algorithm 3: Staggered scheme for MCOP 
initialize ©. initialize @ Va=1,...,N 
initialize boundary condition initialize boundary condition 
t< to t+ to 
while t < fena do while t < feng do 
loop loop 
solve lin. momentum balance solve lin. momentum balance 
// cf. Table 6.1 // cf. Table 6.1 
for a = 1 to N do 
solve evolution eq. for de // cf. solve evolution eq. for d* 
Table 6.1 // cf. Table 6.1 
end 
adapt mesh adapt mesh 
static — AR < Ep static — p maz y [ġe l < Epe 
if static then if static thea 
tet+At tt+At 
increment boundary increment boundary 
condition condition 
break break 
end end 
end end 


end end 


Numerical discretization The proposed models result in a system of par- 
tial differential equations, consisting of the balance of linear momentum and 
multiple evolution equations for the crack order parameters. In this work the 
staggered approach of Miehe et al. [29] is used, which is based on a operator 
split. Each partial differential equation is solved by assuming the other fields 
constant. Together with a time stepping scheme 


ot =o! + Are, ett age" + Arde, (6.26) 


and a time step Af, this results in linear partial differential equations. The 
index n denotes the order parameters for an old time step, while n + 1 denotes 
the order parameter for a new time step. 
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In order to ensure a quasi-static crack propagation, the steady state conditions 
[Pcl < ea; loela < en (6.27) 


are introduced, with the infinity norm |-|,, and a tolerance parameter eg, = 1074. 
After solving the individual equations in each iteration, the condition is eval- 
uated. Only if the condition is fulfilled, the system will progress further in 
time. As a consequence there is an inner iteration loop, which solves each 
equation subsequently, until a steady state is reached, and an outer loop, which 
can be associated with a physical time and a time-dependent load, using bound- 
ary conditions. The procedure corresponding to this approach is illustrated as 
pseudocode in Algorithm 2 for the SCOP model and compared to the MCOP 
model (Algorithm 3). The mobilities M and M“ also reduce to numerical pa- 
rameters and will not influence the results within a certain reasonable numeri- 
cal range. In addition, an exemplary study of the evolution of these staggered 
systems and their iterations is conducted in A.3. In order to solve the partial 
differential equations, a FEM approach, based on the C+ finite element library 
DEAL.II [163], is used*. In this work, first-order Lagrange finite elements and 
a second-order Gauss-Legendre quadrature rule for numerical integration are 
used. The evolution equations are solved using an implicit time stepping pro- 
cedure (cf. Equation (6.26)), except for the elastic driving force term, which is 
considered explicitly in time. Adaptive mesh refinement (AMR) is used to re- 
duce the computational effort. Especially for the setup chosen in section 6.4.2, 
big domains are required which are computationally unfeasible on a uniform 
mesh. In this work a basic AMR strategy is used, and the reader is referred 
to, for example, Heister el al. [150] for a more state-of-the-art AMR approach. 
The strategy of our current approach uses different criteria for refinement or 
coarsening for the evolution equation and the equilibrium of the linear momen- 
tum. Based on the gradient of the order parameter or strain energy density and 
the extent of the elements, discrete changes are calculated for each element. If 


4Results taken from Schiller et al. [1, 3] were conducted based on DEAL.II, and all additional 
phase-field simulation studies used an implementation in the PACE3D framework [164]. 
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certain values are exceeded or undershot, the element is refined, or coarsened, 
respectively. In addition, this procedure is explicitly evaluated in time. De- 
spite limitations of the strategy acknowledged, it enables the computation of 
bigger domains. Combined with a very conservative refinement controlled by 
the values for refinement and coarsening, and frequent execution, the strategy 
has proven its robustness. However, this has the disadvantage of low compu- 
tational optimization compared to more advanced approaches such as those of 
Heister el al. [150]. Nevertheless, it provides huge improvements over a uni- 
form mesh. 


Residual stiffness For a FEM approach, an interpolation function h(6,), 
which will be zero for a fully broken state, will result in a singular stiffness 
and an ill-posed problem. This can be avoided by replacing the degradation 
function by 


3 h(de) 1-0 > os", 
h(d.) = 6.28 
(4) nn else, ) 


with a threshold value oth = 1074. This function will preserve a certain value 


and even ensures a residual stiffness even for completely damaged regions. 


6.3 Extension to a single-obstacle potential 


So far, only a phase-field crack propagation models based on a quadratic po- 
tential term have been introduced. In the context of the classical phase-field 
models, this is often referred to as the single-well potential, or in the context 
of the mechanical community as the AT2 model, cf. e.g [156]. An alternative 
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formulation used the single-obstacle potential, leading to the AT] model. The 
free energy for the SCOP model follows with 


N 
Flu, ġe VO] = f ODIA (6.29) 
+1 y gege &.|v® et dv 
24 cS ice ned | ak 


where the prefactor 9/16 ensures the representation of the corrected surface 
energy of a corresponding sharp crack surface. Analogously, the free energy 
for the MCOP model can be rewritten as 


Fug V9) = | noose (6.30) 


+508 (ea vor + azt) dv. 

Note that the use of a single-obstacle potential leads to a constant term in the 
evolution equation. To ensure that @.,@% € [0,1] Va = 1,...,N, further tech- 
niques such as the active set method or the augmented Lagrangian method, 
cf. e.g [149, 150], have to be used. However, the single-obstacle variant has 
some advantages. Since the derivation of the potential term does not vanish 
for zero damage, the analytical profile has a finite width. Together with the 
well-defined obstacle against initial crack growth, this leads to more localized 
crack growth and avoids damaged materials even far away from the from the 
crack tip. 


In the following work, both the single-well and the single-obstacle variants 
will be used. Since the single-well potential is well established in the liter- 
ature, most of the basic investigations will be performed based on it, if not 
stated otherwise. Subsequently, e.g., some applications to FRPs will also be 
performed with both models or, due to the numerical advantages, only with the 
single-obstacle model. 
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6.4 Results 


Effective crack order parameters® In order to be able to compare the 
results of the SCOP and MCOP model, additional effective order parameters 
are introduced. An order parameter d* can be decomposed into an effectively 
damaged part 6% and an effectively undamaged part $2: 


o% = de + 02. (6.31) 


Each of these parts is defined using the individual crack order parameters of 
the MCOP model 


0S = 97S, oY = 6% (1-62). (6.32) 


Also a totally effective crack order parameter is formulated: 
~ N ~ 
d= YOY. (6.33) 


Note that these effective quantities describe the damaged and the undamaged 
fraction with respect to Er #% = 1, which is in contrast to #, which only 
describes the ratio of damage with respect to 6%. 


6.4.1 Steady-state profiles in 1D® 


To illustrate the difficulties of the heterogeneous SCOP model, steady-state 
profiles are examined for different model parameters and compared to the an- 
alytical solution. For the sake of simplicity, the system is assumed as binary 
and one-dimensional, while no mechanical loads are considered. Instead, a 
completely damaged point, in the middle of the domain, is placed and imposed 


5The content of this section has been taken directly from Schiller et al. [1] with minor linguistic 
changes. 
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x/Eg, 


Figure 6.4: Steady state profiles for a diffuse interface at x = 0, with G?/G! = 100 and Ep = Eps- 


Figure 6.5: (a) Error of the total crack energies eg,. (b) Maximum of the local error of the energy 
density ep (b) and a diffuse interface with varying model parameters. 
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via the boundary condition 6. = 1 or @& = 1, respectively. In addition, an 
interface with different critical energy release rates for both materials is also 
placed in the middle of the domain. In Figure 6.4, exemplary steady-state pro- 
files &., for the SCOP model, and 6, for the MCOP model, are shown along 
with the imposed interface. Note that a high contrast of the crack surface ener- 
gies G2/G! = 100 and the same length parameters Ep = Ep, are chosen. This 
results in a highly distorted profile for the SCOP model, compared to the ana- 
lytical profile, due to the additional spatial gradient in ee o°G*. In contrast, 
the evolution equations of the MCOP model are independent, as the coupling 
is modeled purely by the degradation of the strain energy density, which is not 
present here. Together with the constant G®, this will avoid any distortions, 
and the analytical profile is reproduced. 


The distorted profiles result in different energy densities and do not yield the 
analytically desired crack surface energy G3 = 1/2 (Gl + G2) when integrated 
over the domain. For the investigation of these deviations in the total energy, 
the error estimator deviation 


E_(* 
20 = (6.34) 


eG. = G: 
c 


is introduced. For varying ratios of the critical energy release rates and inter- 
face widths, the deviations are displayed in Figure 6.5a. Due to the higher 
spatial gradient in the critical energy release rate, an increasing deviation for 
a higher contrast G2/G!, can be observed. For Ep, > €g,, the change of the 
critical energy release rate is distributed over a larger physical width, leading 
to lower spatial gradients and lower deviations. If the case &, < Ep, is consid- 
ered, the spatial gradient increases, but its influence is limited to only a part of 
the crack interface, resulting in a relatively low deviation for the total energy 
of the system. From a numerical point of view, £, ~ €, is desired to reduce 
the effort of resolving both interfaces. In this case, however the deviation is 
rather high. The deviations shown here do not exceed ~ 7.8%, which seems 
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tolerable. Nevertheless, the maximum local deviation of the energy density can 


be defined by 
ef= et U = (6.35) 


with the analytical energy density f*, defined by the analytical profile, and 


X=Xmax 


is evaluated at xmax, for which the absolute difference of the density is highest. 
For the various parameters, this local deviations is shown in Figure 6.5b, which 
also shows an increasing deviation with higher G2/G!. With increasing £9, /Edes 
the spatial gradient in YY @%G% becomes more local, which also increase the 
deviations. For the MCOP model, the error quantities eg, and e will vanish, 
as the profiles are identical to the analytical ones and will reproduce the correct 
crack surface energy Gž. 


Although a high contrast in the surface energies only produces moderate devia- 
tions in the total energy, the local deviations of the energy density is an order of 
magnitude higher. So far, no crack propagation behavior has been considered. 
However, it is expected that the inability of the SCOP model to reproduce the 
energy density and total crack surface energy will also influence the evolution 
of the crack, which will be investigated in the subsequent sections. 


6.4.2 Sloped binary interface® 


This section follows the investigations by Henry [55]. Along with a binary 
interface, an initial crack is placed in a two-dimensional rectangular domain, 
with an angle @ between them, cf. Figure 6.6. Note that the real size of the 
domain is large, compared to the interface widths, to avoid any influence from 
the boundaries. 


©The content of this section has been taken directly from Schiller et al. [1] with minor linguistic 
changes. 
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ex 


Figure 6.6: Exemplary structure of an inclined interface, with an angle @ between the binary in- 
terface and the initial crack. A time-dependent stress is applied on the left and right 
boundary. 


In order to reduce complexity and avoid any influence of the chosen homog- 
enization scheme, the elastic material behavior of both regions is assumed to 
be equal and isotropic with A! = A?, u! = u?. In his work, Henry imposed 
displacement boundary conditions, but in order to use the LEFM theory in the 
present work, stress boundary conditions are employed to apply an initial mode 
I crack opening. For a straight crack the stress intensity factor results from 


Ki=AovVa, (6.36) 


with the generally unknown prefactor A, the crack length a, and the stress o, 
far away from the crack, which hence is associated with the stress vector t of 
the boundary condition, cf. Figure 6.6. If the stress is increased monotonously, 
this results in crack propagation above a critical stress, followed by an unstable 
crack growth. To be able to investigate the behavior of the crack propagation 
models, a quasi static crack growth with Kj = Ki. is desirable for a constant 
critical value Kj,. To obtain such a stable crack growth, the applied stress 
will be increased after each time step. If the crack propagates, measured by an 
enlargement of the domain, @, = 1 holds, the stress at the boundaries is reduced 
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Figure 6.7: Simulation domain with an initial crack and a marked and enlarged artificial interface 
region (a). The crack order parameter d. for the SCOP model (b) and ĝe for the MCOP 
model (c), with contour lines for 6 and 62. In addition, the solid boundaries of the 
interface are shown. 


below an estimated new critical value, which is based on the crack growth Aa. 
Furthermore, the stress increment for each time step is also reduced. This 
approach aims to achieve linear crack growth over time and avoid unstable 
crack growth. 


Artificial interface Firstly, the interface between the region is considered 
as artificial. Hence, the crack resistance is assumed to be equal, G! = G2, to- 
gether with the same stiffness, both regions can be considered to consist of the 
same material. The ratio of the interface width parameter of the solid and the 
crack is €9, /l; = 2/5, and for the interface angle, @ = 50° is chosen. The results 
after simulating a mode I crack propagation of such a system are shown in Fig- 
ure 6.7: It could be observed that the initial crack propagated straight through 
the domain, despite the presence of the interface. For the MCOP model, con- 
tour lines of the effective parts ĝl and 2 are displayed, additionally. In the 
solid interface, the crack is transferred from one region to another, continuing 
the total effective crack, consisting of the sum of @%. The resulting effective 
crack Ês of the MCOP model is comparable to d. of the SCOP model. In Fig- 
ure 6.8a, the temporal evolution of the crack length a, normalized by the initial 
crack length ao, is plotted against the normalized simulation time. Here, both 
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Figure 6.8: Normalized temporal evolution of the crack length (a), and the stress at the boundary, 
applied against the crack length (b), for the SCOP and MCOP model with 0 = 50°. 


models also showed the same effective behavior, and the transfer of the crack 


propagation, from one crack order parameter 0% 


to the other, can be observed. 
The evolution of the applied stress o, normalized by the critical stress to prop- 
agate the initial crack Oo, is displayed in Figure 6.8b. Because of the way the 
boundary condition is applied, it should be noted that only the time step is 
shown in which the crack grows, while the other steps are omitted. The mod- 
els show similar profiles and coincide with the expected profile from LEFM, 
cf. Equation (6.36). Regardless of the use of multiple crack order parameters, 
the novel MCOP model is able to reproduce the crack path and kinetics of the 
SCOP model. Furthermore, the presented approach ensures quasi-static crack 
growth and produces an almost linear crack growth over time, as well as the 
desired relation between applied stress and crack length. 
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Figure 6.10: Cross-section with the underlying solid interface, along the binary interface of the 
SCOP and MCOP model, shown in Figure 6.9, for a sharp (a) and diffuse interface 
(b). 


Infinitely tough region The other limiting case of crack propagation with 
a binary interface is an infinitely tough interface, which can be obtained when 
G2 = œ, which ensures that the upper region cannot be damaged at all and 
enforces the crack to grow along the interface. For both regions, the same stiff- 
ness parameters are used, while the boundary condition is applied as before. 


In Figure 6.9, the crack path for both heterogeneous models and @ = 50° are 
shown, considering both a sharp solid interface (a-c) and a diffuse solid inter- 
face (d-f). The marked cross section is also drawn in Figure 6.10 along with 
the solid interface. The SCOP model shows a non-physical behavior in some 
aspects: 
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e Due to the interpolation of the individual critical energy release rates G% 
and the infinitely tough upper region, no crack can occur there or in the 
interface region. This constrains the crack propagation in such a way 
that the crack is deflected before it reaches the solid interface. 


e This interpolation also causes problems as to how the crack grows along 
the interface. Because the order parameter in the interface must be zero, 
this leads to a certain distance of the crack from the interface. 


e For the diffuse interface, the high gradients of @, towards the interface 
combined with the nonconforming mesh cause numerical issues, lead- 
ing to numerical artifacts, as shown in the magnification of Figure 6.9b. 
Regarding the sharp interface, which requires a conforming mesh, these 
artifacts are not present. 


In contrast, the MCOP model exhibits fewer of these difficulties: For the dif- 
fuse interface, the crack grows straight up to the interface, where it performs 
distinct deflection. The deflection is less pronounced for the sharp interface 
due to the loss of the driving force immediately beyond the interface. However, 
for both the diffuse and sharp interfaces, the crack continues directly at the 
edge of the interface, and the crack exhibits the desired analytical profile of the 
diffuse interface, as illustrated in Figure 6.10. 


For a quantitative comparison of the models, the approach of Henry [55] and 
the analytical analysis of Amestoy and Leblond [13] are used and briefly intro- 
duced, in the following. The stress intensity factors for the crack modes I and 
II of the straight crack, i.e., before the deflection of the crack at the interface, 
results from 


Kı = Aooya, Ky = 0. (6.37) 
In contrast, the stress intensity factors can be described by 


Ki = f ()Kio, Ku = 8()Kio, (6.38) 
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Figure 6.11: Exemplary visualization of the evaluation of the crack propagation (a) and compar- 
ison with an analytical solution for the SCOP and MCOP model (b), for a crack 
propagating along an infinitely tough interface, with a different angle 0. Various 
solid length parameters &,, as well as a sharp interface (SI) are examined. 


directly after the kink of the crack path, where f(@) and 8(@) are given by 
Amestoy and Leblond [13]. Ajo is the stress intensity factor right before the 
kink of the mode I crack. The energy release rates after the kink G results from 


22, £2 
Treen 


Be (6.39) 


while the ratio of G and the energy release rates before the kink Go can be 
expressed using the analytical solution 


G 
u Ce Ce (6.40) 
0 
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with E’ = E /(1 — v?). With Equation (6.37), this ratio can also be given by 


ô (ô\ 

e (=) : (6.41) 
using the applied stress at the interface after the kink 6, and a reference stress 
Oo, representing the stress of a straight crack with the same crack length in the 
y-direction. In Figure 6.1 1a, an exemplary illustration of the evaluation proce- 
dure is given. When the crack hits the interface, the critical stress for propaga- 
tion increases, which can be observed by a peak in the graph. By determining 
the maximum value of this peak for which a, > ayọ holds, and comparing it 
with the reference stress 6 of a simulation without a solid interface but with 
the same crack length in the y-direction a,, a comparison with the analytical 
solution is possible. In Figure 6.11b, this procedure was conducted for both 
models, several angles, sharp and diffuse interfaces. For small angles, both ap- 
proaches show a fairly good agreement with the analytical solution. For higher 
angles, however, the SCOP model strongly differs from the solution, while the 
MCOP model still agrees with the analytical solution from the LEFM. Note 
that a variation of the width of the solid interface was also conducted, where 
the sharp interface can be considered as a limiting case. Neither increasing nor 
decreasing of the interface width has any significant effect on the results. For 
a mesh convergence study of this setup the reader is referred to Section A.3. 


In conclusion, regarding the SCOP model, neither the crack grows to the in- 
terface, nor is the model able to reproduce the stress increase at which crack 
propagation occurs along the inclined interface. In contrast, the MCOP model 
propagates to the interface, develops more significant deflection, and predicts 
the increase in stress according to theory. Regarding the comparison with 
sharp solid interfaces, no significant differences could be found. This limits the 
novel MCOP model not only to the application at diffuse interfaces, but also 
allows for the application to sharp interfaces. Concerning the diffuse bound- 
ary approach, both SCOP and MCOP avoid the need for conforming meshes. 
However, the MCOP model offers the possibility to extend the homogenization 
scheme, e.g., by considering jump conditions [22, 90]. 
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Figure 6.12: Schematic of a single inclusion problem. A square plate with side length Z, a spheri- 
cal glass fiber (GF) inclusion with radius R, is embedded in a thermoset (TS) matrix. 
In some studies, an initial crack is placed at the boundary. 


6.4.3 Single inclusion problem 


To motivate the use of mechanical jump conditions, cf. Section 2.3.3, a sin- 
gle inclusion problem (SIP) is considered. Therefore, a square plate with a 
spherical glass fiber (GF) inclusion is embedded in a thermoset (TS) matrix, 
cf. Figure 6.12. The fiber is represented in the context of a phase-field method, 
i.e. a diffuse interface describes a smooth transition from the matrix to the in- 
clusion, cf. Section 2.3. The strain energy densities of such a system under 
different loads and different homogenization schemes are then investigated. Fi- 
nally, an crack propagation simulation with an initial crack is performed with 
the single-obstacle and single-well and both homogenization schemes. 


6.4.3.1 Strain energy densities 


First, the effective strain energy density f., and the phase-specific strain en- 
ergy densities of both phases are discussed. Therefore, the SIP without initial 
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Figure 6.13: The effective and the phase-specific strain energy densities of a single inclusion prob- 
lem with displacement boundary conditions for the Voigt-Taylor scheme (a) and con- 
sidering jump conditions (b). The diffuse interface is highlighted in gray. 


cracks is considered, cf. Figure 6.12. All following evaluations of the quanti- 
ties are performed from the center of the inclusion along the diagonal of the 
square domain. The elastic material parameters for both materials are listed in 


Table 6.2. 


Displacement boundary conditions In the two-dimensional domain, a 

plane strain state is assumed. In addition, the orthogonal displacement com- 
ponents at the boundary are constrained to yield 1 % macroscopic strain. Fig- 
ure 6.13 shows the strain energy densities for the Voigt-Taylor scheme and 

when considering jump conditions. In both cases, the effective value is lower 
in the inclusion due to the higher stiffness of the inclusion, and a smooth dif- 
fusion to the higher value of the matrix can be observed. When considering 
jump conditions, the phase-specific strain energy density of the inclusion is 
lower than the effective value, as well as the value of the matrix, throughout 
the diffuse whole interface. In contrast, the Voigt-Taylor scheme leads to a 
large increase in the strain energy density of the inclusion towards the outer 
boundary of the diffuse interface. The values of the matrix is also in the or- 
der of magnitude lower than that of the inclusion. Since the phase-specific 
strain energy densities are part of the driving force for crack propagation in the 
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Figure 6.14: The total effective and the phase-specific strain energy densities of a single inclusion 
problem with eigenstrains in the matrix for the Voigt-Taylor scheme (a) and consid- 


ering jump conditions (b). The diffuse interface is highlighted in gray. 


context of the MCOP model, the use of the jump condition is preferable. In 
this case, for example, assuming the same crack resistance, crack nucleation 
of the inclusion would be favorable for the Voigt-Taylor scheme. For the jump 
condition a crack nucleation of the matrix occurs only for a significant higher 


load. 


Eigenstrains in the matrix In addition to the external load due to a dis- 
placement boundary condition, the system is extended to include an eigenstrain 
in the matrix material, e.g., to mimic curing shrinkage. Therefore, the orthog- 
onal displacement components at the boundary are set to zero and an inelastic 
strain of —2 % of the matrix is assumed. As before, Figure 6.13 shows the 
strain energy densities for the Voigt-Taylor scheme and if jump conditions are 
considered. In both cases, the qualitative behavior is similar. Towards the outer 
boundary of the interface, the energy density of the inclusion exceeds that of the 
matrix. Also, the highest effective value is at a similar distance from the center. 
A quantitative comparison reveals that in particular the strain energy density 
of the glass fiber is predicted to be more than 20 times higher compared to the 
jump condition framework. In the context of the MCOP model, this would lead 


to vastly different eigenstrains required to initiate crack growth. 
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6.4.3.2 Crack propagation 


Next, the setup is extended to crack growth, so an initial crack is placed in the 
domain, cf. Figure 6.12. In addition, the eigenstrain, which mimics the shrink- 
age during curing, is applied linearly in time, allowing the single-obstacle and 
single-well potentials, as well as the Voigt-Taylor and jump condition schemes, 
to be compared. The elastic material parameters for both materials are listed in 
Table 6.2 and the crack resistances are assumed with GIS = G® = 50Jm~?. 


Boundary and initial conditions A decomposition of the strain of the 
matrix, the thermoset (TS), is defined by 


eS =ef +e, (6.42) 


cf. Section 5. Therefore, ey accounts for a volume change due to curing and 


is prescribed in this section by 


(6.43) 


where x"° is initially zero (TS (t =0) = 0) and decreases linearly with time to 
mimic curing shrinkage. In addition, an effective plane strain state and periodic 
boundaries are assumed and thermal effects are neglected. 


Results Figure 6.15 shows the effective maximum principal stress, cf. Sec- 
tion 5.21 and Equation (5), versus the shrinkage of the thermoset. The evo- 
lution is depicted for the single-well and single-obstacle potentials, as well 
as for the Voigt-Taylor scheme and considering mechanical jump conditions 
in both cases. The single-well potential shows the typical stress degradation 
associated with such a model. Even without propagation of the completely 
broken state, a strong degradation of the effective behavior occurs. In contrast, 
a sharp decrease in stress is observed for the single-obstacle potential. Unsta- 
ble crack growth is characterized by an almost stepwise decrease in the stress 
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Figure 6.15: Effective maximum principal stress versus the curing shrinkage of the thermoset, 
cf. Equation 6.43 for the single-well and single-obstacle potential, using the Voigt- 
Taylor scheme and considering mechanical jump conditions. In addition, individual 
time steps are annotated with numbers and displayed in Figure 6.16. 


of the corresponding curve and some time steps are labeled with individual 
numbers. For the latter, the effective crack order parameters are shown in Fig- 
ure 6.16, although the discussion is restricted to the single-obstacle simulation 
cases. After the initial crack has grown straight, the interfacial crack stops at 
the stiffer inclusion, which is an obstacle to direct crack propagation (1,2). In 
both cases, the Voigt-Taylor scheme and the mechanical jump condition, the 
crack grows around the inclusion and merge (3,4). However, the Voigt-Taylor 
scheme exhibits crack growth inside the inclusion, leading to a more triangular 
crack path, while considering the jump condition leads to a round path outside 
the inclusion. Subsequently, with the Voigt-Taylor scheme, the cracks grow 
perpendicular to the initial crack along the cross section with the least matrix, 
crossing the periodic boundary and finally stopping again at the inclusion (5). 
For the jump conditions, the same curing shrinkage x"° is not sufficient to ob- 
serve any crack growth perpendicular to the initial crack beyond the merging 
of the crack tips (6). 


Figure 6.17 summarizes the final effective crack path ĝe for 77° = —1.75%, as 
well as the crack order parameters for both materials #75, °F. As discussed 
previously, the use of a single-well potential also leads to an increase in the 
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Voigt-Taylor 


Jump conditions 


Figure 6.16: Effective crack path for single-obstacle potential, using the Voigt-Taylor scheme, and 
considering mechanical jump conditions at different time steps, see Figure 6.15. The 
sharp representation of the inclusion is shown as white contour lines. 


crack order parameter far away from the fully broken crack, resulting in a re- 
duced effective stiffness, cf. Figure 6.15. Apart from this characteristic, the 
final crack paths are similar. Much more noticeable is the influence of the ho- 
mogenization scheme used, since in addition to Figure 6.16, the contributions 
to the effective crack order parameter can also be discussed. In all cases where 
the mechanical jump condition is considered, only the thermoset is damaged, 
resulting in a circular crack path around the inclusion. The use of the Voigt- 
Taylor scheme results in a partially completely damaged interface. Since the 
crack resistances are the same but the Young’s moduli differ substantially, such 
behavior is unexpected in the context of the K concept, cf. Section 2.4. Finally, 
the crack path is growing perpendicular to the initial crack in the Voigt-Taylor 
cases, while the magnitude of the curing shrinkage is not high enough for such 
a change in crack growth for the jump conditions cases. 
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Figure 6.17: Final crack paths for single-obstacle potential, using the Voigt-Taylor scheme and 
considering mechanical jump conditions. The effective order parameter d. as well 
the phase-specific variants for thermoset @7S and glass fiber 6" is displayed. 
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Elasticity Fracture 
Material E v Reference Ge Parameter Value 
Thermoset 3.45GPa 0.38 [165] 100.0 Jm~? l, lum 
Glass fiber 73.0GPa 0.22 [166] 200.0 Jm~? Ep 1.25 um 
(a) (b) 


Table 6.2: Material properties (a) and interface width parameters (b) for simulating the fracture of 
FRP volume elements. 


6.4.4 Application to FRPs’ 


In Section 6.4.2, the same elastic material behavior has been assumed for both 
regions. But realistic systems with a high contrast in crack resistance will 
most likely also have a high contrast in their elastic parameters. In this section, 
a fiber-reinforced polymer (FRP) is chosen to schematically demonstrate the 
ability of the novel MCOP model to illustrate crack propagation behavior in 
the context of a material with hetergeneous elastic properties. A quantitative 
analysis of the results is omitted, as this would most likely require an extension 
of the model, e.g., a state-of-the-art tension-compression splitting, or account- 
ing for interfacial crack propagation, which is beyond the scope of this work. 


The matrix consists of a UPPH resin system [167], reinforced by glass fibers. 
The material parameters for both materials are given in Table 6.2a, while the in- 
terfacial widths used are shown in Table 6.2b. FRPs exhibit a complex fracture 
behavior: Either the matrix may fail, the fibers may break, or the material may 
fail due to the debonding of the interface. To investigate crack propagation in 
such a material, volume elements with a certain fiber volume content, orien- 
tation, and periodicity are considered. Boundary conditions, such as normal 
Neumann or Dirichlet types, do not account for the periodicity of the domain. 


7The content of this section has been taken directly from Schöller et al. [1, 3] with minor linguistic 
changes. 
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Hence a periodic type is chosen: In addition to the periodic order parameters 
and displacement fields, a superimposed periodic displacement boundary con- 
dition [168] is applied in such a way that the macroscopic strain of the volume 
element follows by 

€=e£(t)e,®e,, (6.44) 


where the normal strain in the x-direction, €(t), will be increased linearly with 
time, until the volume element exhibits a complete failure. 


Unidirectional 2D volume elements Firstly, unidirectional fiber rein- 
forced volume elements are be investigated. This makes it possible to reduce 
the system to a 2D system, so as to reduce the computational effort. The square 
volume elements with a side length of 100m contains fibers with a radius of 
4um. For different fiber volume fractions vy, various realizations will be con- 
sidered in the following. In Figure 6.18, the effective crack d. and the contour 
lines of the fibers are presented in a fully broken state. All realizations show 
an overall crack direction, perpendicular to the applied load. Thus, the desired 
and dominant crack opening mode I is reproduced, and the crack paths tend to 
become more complex, when using a higher fiber volume fraction, as the fibers 
become an obstacle for the crack, which results in an elongation and contortion 
of the crack path. This is primarily caused by the lower crack resistance of the 
matrix. Thus, the crack also often propagates through matrix-dominated re- 
gions, but also between fibers caused by stress concentration. Since the MCOP 
model does not account for the failure mechanism of interfacial debonding, 
this failure mechanism cannot be observed when the crack propagates in fiber 
dominated regions. Due to the periodic boundaries, the complete failure of the 
volume element forces the crack tips to merge. After the merge, partly ‘dead 
ends’ can thus be observed in some realizations presented in this work. 
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vf Different realizations 


Figure 6.18: Different realizations of randomly generated periodic 2D volume elements with uni- 
directional fibers for varying fiber volume fractions vf. Additionally, the effective 
crack &., after uniaxial strain and failure is shown. 
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(a) (b) 


Figure 6.19: A randomly generated periodic 3D volume element with isotropic fiber distribution 
(a) and the fractured volume element after applying an uniaxial strain, highlighted by 
ared crack surface (b). 


Isotropic 3D volume element In the case of a 3D volume element, the 
fiber volume content vf = 25 % is chosen with an isotropic orientation distribu- 
tion. The latter implies that there is no preferential direction in the fiber distri- 
bution. As before, the volume element has a side length of 100 um, consisting 
of fibers with a radius of 4um and a length of 80um. The volume element 
was generated using the approach of Schneider et al. [133] and is shown in Fig- 
ure 6.19a. As in the 2D case, a macroscopic strain is applied in one direction, 
in addition to the periodicity, and is increased linearly over time. The failed 
domain with the red crack surface is shown in Figure 6.19b. As for the 2D 
simulations, the crack surface is mainly perpendicular to the load direction and 
occurs solely in the matrix, but still shows a quite complex crack path, due to 
the fiber distribution and the stress concentration that arise from it. 


Considering mechanical jump conditions In this section, a glass fiber- 
reinforced polymer is chosen to schematically demonstrate the advantage of 
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Figure 6.21: Detail of the crack path from Figure 6.20f (Voigt-Taylor, GEF /GTS = 1/2) for the 
effective crack order parameter ĝe (a), the thermoset crack order parameter or S (b), 
and the glass fiber PC (c). 


considering mechanical jump condition compared to the application of the 
Voigt-Taylor scheme in the context of MCOP. This material system is moti- 
vated by the high contrast regarding the elastic properties of the matrix mate- 
rial, a thermoset, and the glass fiber: A Young’s modulus of ETS = 3.45 GPa, 
and ECF = 73.0GPa is assumed for the glass fiber, respectively the thermoset. 
The Poisson’s ratio of fiber and thermoset are vTS = 0.38 and v® = 0.22 ac- 
cording to [165, 166]. For simplicity, unidirectional fiber-reinforced volume 
elements are used. This allows a reduction to a 2D system. A square with a 
side length of 100 um and a volume fraction of 40 %, with a fiber radius of 
4 um is chosen. As boundary conditions, the macroscopic strain tensor 


E(t) =Eu(t)e,®e,, (6.45) 


is applied. The function &,.(r) will be increased linearly with time, until the 
volume element fails completely. While the crack resistance of the thermoset is 
kept constant with GTS = 100.0Jm~*, the crack resistance of the glass fiber is 
varied. Despite a failure of the fiber being possible under certain circumstances, 
a failure of the matrix is expected for all ratios used in this work. With the 


133 


6 Phase-field modeling of crack propagation in heterogeneous materials 


assumption that the crack surface energy will be minimized and GTS > GO, 
failure of the fiber seems preferential. But even for such a case, the matrix can 
exhibit a lower crack resistance in the sense of stress intensity factors. So not 
only the crack surface energy, but also the capability to provide this energy, 
e.g., by a higher Young’s modulus, determines which material component fails. 
In Figure 6.20 the final crack paths, with the effective crack order parameter 
field 

6. = 675 TS | gGF GF (6.46) 


are displayed. Regarding the simulations that account for the jump conditions 
(a-c), the system exhibits the same crack path for all variants. From a nucle- 
ation of the crack between close fibers due to stress concentration, a matrix- 
dominated path is predicted. For the Voigt-Taylor scheme, the two higher ra- 
tios (d,e) show the same crack path, even if more fiber damage is present. For 
the lowest ratio (f) a completely different path is predicted, with damaged fiber, 
even far away from the final failure crack path. In addition, the crack propa- 
gates at the inner side of the fiber interfaces, which is in contrast to the expected 
behavior. In Figure 6.21 a detailed section of the Figure 6.20f is displayed. In 
addition to the effective crack order parameter field (a), also the crack order 
parameters fields of the thermoset (b) and glass fiber (b) are provided. In the 
latter, it can be observed that in addition to the matrix material also the fiber 
exhibit completely failure, which is in contrast to the expected behavior. 


The macroscopic stress-strain curves are displayed in Figure 6.22. Already in 
the linear regime of the curves, a difference can be observed: As the Voigt- 
Taylor scheme describes an upper limit for the elastic energy, it seems reason- 
able that it exhibits higher stresses. Since it is only accurate for specific cases, 
such as a parallel material chain, the difference in stress is in accordance to 
theory [60]. As the failure is dominated by the matrix behavior, the influence 
of the varied crack resistance of the glass fiber is negligible. However, a small 
influence during the crack nucleation between two fibers due to the diffuse 
interface remains. On the contrary, the failures for the Voigt-Taylor scheme 
occurs at quite different strains. The assumption of same strains for fiber and 
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Figure 6.22: Macroscopic stress-strain curves of 2D unidirectional fiber-reinforced volume ele- 


ments for a Voigt-Taylor scheme and if jump condition are considered, with (y) = 
+ Jo wdv. In addition, the ratio of the crack resistances GCF /GTS is varied. 


matrix leads to high driving forces and therefore an earlier failure of fiber and 
matrix, compared to the case where jump conditions are considered. 


6.5 Interim conclusion 


In this chapter, different phase-field models for crack propagation in heteroge- 
neous systems have been presented. The SCOP model, which uses a single 
crack order parameter to account for damage and is based on established ap- 
proaches, and the novel MCOP model, which uses multiple order parameters 
to account for damage in each region separately. In combination with the jump 
condition approach of e.g. Schneider et al. [61], the novel model is able to 
provide a more accurate representation of the driving forces for crack propa- 
gation in heterogeneous materials compared to established crack propagation 
phase-field models. It has been shown that the MCOP model is able to repro- 
duce the surface energy of the sharp interface for an interfacial crack, while 
the SCOP has a high error in the local energy density distribution over the 
diffuse interface. In addition, a study of crack propagation along an sloped 
interface indicates a higher qualitative and quantitative accuracy of predicting 
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the crack growth than established models when compared to an analytical so- 
lution. Based on a single inclusion problem, the advantages of considering the 
mechanical jump condition in the MCOP were demonstrated and the single- 
well and single-obstacle potentials were compared. Moreover, the application 
of the MCOP model to FRTS was subsequently shown, demonstrating the ad- 
vantages of the model for such systems, but not limited to this material class. 
These simulations primarily use boundary conditions to apply the loading, such 
as a uniaxial tensile test. In contrast, in the following chapter, the curing model 
from Chapter 5 will induce inelastic strains that cause damage to the FRTS 
during curing. 
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In this chapter, the curing simulations of Chapter 5 are combined with the 
MCOP crack propagation model of Chapter 6 to predict crack formation dur- 
ing the manufacturing process. Therefore, similar to previous simulation stud- 
ies, crack formation in different levels of abstraction of real microstructures of 
FRTS are considered and compared. In addition, a modified approach for the 
choice of an optimal mobility to reduce the computational effort is introduced. 


7.1 Optimal choice of mobility 


The introduced phase-field crack propagation models are solved on the basis 
of an Allen-Chan type equation, cf. Section 6.1.1, which yields evolution equa- 
tions for each crack order parameter. Together with the staggered approach, 
these evolution equations and the linear momentum balance are solved. For 
the implementation in the context of DEAL.II [163] , these PDEs are solved us- 
ing a time implicit FEM approach. In contrast, the PACE3D frameworks [164] 
relies on an explicit time integration scheme for the order parameters. In the 
context of the multiphase-field method, this allows fast execution even for a 
large number of phases using e.g. local reduced order parameter techniques. A 
disadvantage of the scheme is the limited numerical stability of the time inte- 
gration. For crack propagation models, a fixed but otherwise arbitrary mobility 
is usually chosen and a time step width Ar within the stability region is used. 
This approach ensures a stable time integration, but can lead to a large number 
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of inner iterations necessary to satisfy the static crack propagation criterion, cf. 
Section 6.2. Therefore, an optimal choice of mobility is desirable to reduce the 
computational cost, and will be derived in the following. 


A generic evolution equation for a crack order parameter 6° can be written by 
of =AAGS + BOE +C, (7.1) 


with a prefactors A for the diffuse term, B for the linear term and an independent 

term C. In this section the investigation are restricted to the MCOP model with 

the single-obstacle potential, whereby these prefactors can be identified with 
Me 


2M% 9 
A=M°GL, ae C= = (218 De, ct). (7.2) 


ic 


For the implementation in PACE3D a finite volume approach is chosen. To- 
gether with an equidistant grid this leads to a central difference scheme for the 
second derivatives with 


> bi —%i+Gi-1 
oxox Ax2 


(7.3) 


for a derivative in the spatial direction x with a cell size Ax in this direction. 
Based on this, the CFL condition [169, 170] 


3 
Ax; 
NSS D Peg 


i 


(7.4) 


can be formulated, with i = 1,2,3 for the corresponding direction. With Equa- 
tion (7.2), and accounting for all spatial direction 


3 
AX; 
M Ar<), TEN Re (7.5) 
i 2( E ae) 


follows. This condition is usually used to determine a stable time step width 
At. Instead, the time step width can be kept constant, to e.g. enable a stable 
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coupling of other equations such the heat equation, and a stable mobility M* 
could be chosen. 


However, instead of choosing a global value for the mobility in this chapter, the 
value is determined locally for each cell in the grid. Ultimately, this approach 
results in spatially varying mobilities M“ = M* (x), which is are longer a nu- 
merical parameters to chose freely, but are determined by the equality case of 
Equation (7.5) with a certain safety factor. Note that the finite volume method 
used usually preserves fluxes. With varying mobilities this characteristic is 
no longer valid, but since this is not a necessary restriction for the numerical 
scheme in this work, the chosen approach can be considered as legitimate. 


7.2 Results 


Material parameters The elasticity, thermal, curing and fracture material 
parameters for the UPPH resin and glass fiber are chosen according to the pre- 
vious simulation studies, cf. Chapter 5. Therefore, most of the resin properties 
are taken from Schwab and Denniston [76]. A complete list of references for 
each property is given in Table 5.2. In addition, the material parameters for 
the glass fiber are listed in Table 5.1. Additionally to these already introduced 
material parameters, crack resistances G“ are also required and are assumed to 
be the same as in the previous crack propagation studies in FRPs and are listed 
in Table 6.2a. 


Boundary and initial conditions For the boundary and initial conditions 
for the displacement and temperature field, the same approach is used as in 
Chapter 5. Therefore, the compression molding process is represented by two 
steps, a heating from room temperature to 373 K for 3 minutes and a subsequent 
cooling to room temperature, whereas the curing process takes place during 
the first phase. The volume elements are considered to be periodic, except 


139 


7 Crack formation during curing of FRTS 


for one direction which is associated with the molds through special boundary 
conditions, cf. Section 5.2. 


For the crack order parameters, an initial crack is placed to avoid cumbersome 
techniques to consider crack nucleation in a quantitative manner, cf. Sec- 
tion 6.2. Therefore, a crack with 25 um length is placed in a resin-rich re- 
gion. Although this approach lacks sophistication, it is considered valid in the 
context of this work. The results presented are based on generated volume 
elements and a general statement is possible only based on many of such en- 
sembles. Although several samples are used in this work, they are not enough 
to be statistically significant. Nevertheless could a manual placed initial crack 
still give insights to the basic fracture behavior of FRTS. Moreover, such an 
initial crack can also be associated with imperfections such as voids or a not 
completely closed flow front resulting from an uneven flow of the initial charge 
through the geometry of the final component. In this sense a manual placement 
of an initial crack seems admissible. 


Unidirectional fiber-reinforced thermoset First, unidirectional fiber- 
reinforced volume elements in a 2D system are considered. Figure 7.1 shows 
the final crack path for 10 %, 30% and 50% fiber volume content for three 
different samples. In addition, the initial crack is shown in red contour lines 
and the corresponding sharp fiber-resin interface is displayed in white contour 
lines. 


Although there is no preferred crack growth direction in the microstructure, the 
direction of the initial crack also determines the preferred direction of crack 
growth during cure. Thus, a straight final crack from the initial crack is ex- 
pected for a pure resin system. For 10% fiber volume content, only a small 
deviation from such a straight crack path can be observed. This occurs because 
the fiber represents an obstacle to straight crack propagation. Some ensemble 
members, such as Figure 7.1d, show a complete failure of the volume element 
in one direction, where the fibers only slightly alter the straight crack path. For 
30 % fiber volume content, the crack paths are more deflected by the fibers. 


140 


7.2 Results 


10% 30% 50% 


Figure 7.1: Final effective crack order parameter field after curing for different fiber volume frac- 
tions and ensembles for the unidirectional fiber-reinforced volume elements. In ad- 
dition, contour lines represent the corresponding sharp fiber-matrix interface and red 
contour lines represent the initial crack. 
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For example, Figure 7.1e shows a strong deflection from the initial crack ori- 
entation during crack growth due to the fibers. Overall, the final crack path is 
less spread over the volume elements and no complete failure occurs. For the 
50 % ensembles, there is no crack growth other than the initial growth. This is 
mainly due to the stiffer fibers, which exhibit no chemical shrinkage and lower 
thermal expansion. In higher amounts, they strongly reinforce the system and 
transfer the main load so that the initial crack is almost stress free. This is also 
illustrated by the maximum principal stress shown on the left in Figure 7.2. 
There exist a circular peak in stress in the fibers surrounding the initial crack, 
while the crack itself has almost no stress and therefore lacks the necessary 
driving force to propagate the crack further. 


Note that the unidirectional fiber-reinforced sample does not accurately depict 
the real microstructure of a compound molded glass fiber-reinforced UPPH 
system. Thus, the fiber cross sections presents a relatively low obstacle to crack 
growth and the effective crack behavior seems to be mainly dominated by the 
volume fraction. Therefore, an extension to a more realistic microstructure 
seems reasonable and is carried out in the following. 


Two-dimensional curved fiber-reinforced thermoset In the next step, 
the microstructure is represented by curved long fibers. In order to reduce 
the computational effort, the simplification to 2D is maintained. As before, 
Figure 7.3 shows the final crack path for 10%, 30 % and 50% fiber volume 
content for three different samples, along with the initial crack and the corre- 
sponding sharp fiber-resin interface illustrated by red and white contour lines, 
respectively. In contrast to the unidirectional fiber-reinforced elements, none 
of the 10 % ensemble members exhibits a crack across the entire domain, nor 
are they straight. Instead, the cracks follow the curved fibers either directly at 
the interface or a deflections before reaching the interface. The latter is likely 
due to the complex stress state and different eigenstrains of the materials. As 
the volume fraction of the fiber increases, the crack path becomes more con- 
strained, resulting in shorter crack paths and more pronounced deflections and 
kinks in the crack path. As before, the 50 % ensemble shows no crack growth 
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Figure 7.2: Distribution of maximum principal stress o* for volume elements with a fiber vol- 
ume fraction of 50% for a unidirectional fiber-reinforced system (left) and a two- 
dimensional curved fiber-reinforced system (right) after curing simulation with a crack 
(red). 


other than the initial growth as the high amount of stiffer fibers transfers the 
load. This is again evident from the maximum principal stress displayed on the 
right side of Figure 7.2, where the peaks are in the fibers surrounding the initial 
crack, which experience only low stresses. 


Three-dimensional curved fiber-reinforced thermoset Figure 7.4 
shows a volume element with 10% fiber volume content. In addition to the 
curved fibers, the crack is shown in color, with the different colors representing 
the crack growth at different time steps, starting from the initial disc-shaped 
crack. Due to the volumetric nature of the crack and the three-dimensional 
strain state, a direct comparison with results from 2D simulations is problem- 
atic. The final size of the crack appears small compared to the 2D results. This 
may be due to the smaller domain size, which tends to restrict crack growth 
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10% 30 % 50 % 


Figure 7.3: Final effective crack order parameter field after curing for different fiber volume frac- 
tions and ensembles for the curved fiber-reinforced volume elements in 2D. In addition, 
contour lines represent the corresponding sharp fiber-matrix interface and red contour 
lines represent the initial crack. 
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Figure 7.4: Volume element with 10% fiber volume content and the final crack in color. The 
different colors depicts the crack growth from the initial state to the finial crack. 


due to the boundary conditions. Nevertheless, the volume elements are more 
representative of areal microstructure of the FRTS. Moreover, a complete fail- 
ure or severely damaged microstructure is not expected since the composite 
molding process is capable of producing fully functional components of glass 
fiber reinforced UPPH resin. Therefore, the three-dimensional result supports 
the assumption that microcrack formation during the curing process is limited 
to local areas of high stress concatenation and does not propagate through a 
larger area of the component. 
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Figure 7.5: The crack length factor, the integral of the effective crack order parameter of individual 
simulations normalized by the initial crack in pure resin, against the fiber volume 
fraction. Distinguished by the level of abstraction of the microstructure (case), and 
confidence intervals are provided. 


Comparison To compare the previous results, the integral of the effective 
crack order parameter of individual simulations, normalized by the initial crack 
in pure resin is plotted against the fiber volume fraction in Figure 7.5. In ad- 
dition, two more volume fractions are considered, and confidence intervals are 
given. Note that a higher crack length factor indicates a longer crack, but does 
not specify the path ofthe crack, i.e., whether the crack is straight or more com- 
plex, e.g., deflected by fibers. For the for 10 % unidirectional fibers, the crack 
length is high, due to the mainly straight cracks. While with higher fiber vol- 
ume fraction, the crack length decreases as the cracks do not propagate through 
the whole domain. For the curved fibers, the crack length for 20 % volume frac- 
tion is even increased because the curved fiber allows more crack propagation, 
as the path is deflected along the fibers resulting in longer crack paths. This 
behavior can be observed up to 40 % volume fraction. While these curved fiber 
volume elements still exhibit significant cracks, the crack length decreases dras- 
tically for the relatively evenly spaced unidirectional fiber volume elements. 


146 


7.2 Results 


IWM-018941 10.0kV 7.0mm x65 BSECOMP 


Figure 7.6: Cross-sectional SEM image capturing the fracture within unidirectional glass fiber- 
reinforced thermoplastic. Reprinted from [171]. 


Note that the reduction to 2D is expected to introduce a systematic error. To- 
gether with other sources of error such as material parameters, the quantitative 
aspect of these results should be treated with caution. However, the qualitative 
behavior and the underlying crack formation mechanism are still considered to 
be representative. 


In Figure 7.6, an SEM image of a fractured cross section of a unidirectional 
glass-reinforced thermoplastic is shown. In the absence of directly compa- 
rable experimental results in literature, a comparison based on this image is 
discussed in the following. Note that a different matrix material as well as a 
different loading, i.e. boundary conditions, are considered. Nevertheless, a 
qualitative comparison of the crack mechanism is possible. The unidirectional 
fiber-reinforced simulations shows a matrix dominated failure. In contrast, the 
SEM image shows a pronounced fiber-matrix interface failure. The crack prop- 
agation model used in the present work does not include a special mechanism 
to account for interfacial crack propagation. Thus, it does not take into ac- 
count the surface energy in such interfaces, which would be released for crack 
growth at these interfaces and affect the Griffith fracture criterion [30]. Such a 
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mechanism could lead to a different crack behavior, as a crack at the interface 
could be favorable. Also the initial crack is placed in the matrix due to the 
feasibility of this approach. In real volume elements, a crack is expected to 
nucleate at the fiber-matrix interface. However, after the initial straight crack 
in the resin, the crack propagates along the fiber-matrix interfaces of the curved 
fibers. This indicated that for larger systems with the possible longer cracks, 
the crack could propagate along the interface similar to the experimental image 
taken from an already highly damaged specimen. A direct comparison of of 
this behavior with experimental results is not possible due to the lack of images 
perpendicular to the fiber cross section. Despite a deficit in the comparability 
of simulations and experiments, a similar crack propagation could be observed. 
Nevertheless, an extension of the model used to include crack nucleation and 
interfacial failure could improve the prediction of crack formation during the 
curing process. 
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This work deals with the modeling of curing and crack propagation during the 
compression molding process of FRTS on the nano and micro length scale. In 
the following, the results of the different aspects of this work are summarized 
and final conclusions are drawn. In addition possible extension of the models 
and overall approach is presented. 


Two-stage polymerization using molecular dynamics! Based on 
molecular dynamics simulations, a simplified but industrially used resin sys- 
tem was extended by adding a fiber surface and sizing layers. The approach of 
Schwab and Denniston [76] was used to model the two-stage polymerization 
of the UPPH resin. In addition, y-MPS was chosen for the coupling agent, and 
further assumptions were made about the structure of the sizing. A systematic 
procedure for the development of a final cured system was presented. Based 
on this approach, evaluations of average quantities during the reactions were 
performed. Moreover, the system was also evaluated along the normal of the 
fiber surface, which provides a spatial analysis of the fiber-sizing-resin system. 


Based on the established REACTER framework, coupling agent monomers un- 
dergo a condensation reaction yielding a distribution of monomers, dimer and 
higher oligomers. Due to a lack of information in the literature, validation of 
this distribution did not seem feasible. Nevertheless, this offers an alternative 
approach to an arbitrarily prescribed distribution. In both the UP as well the 


'The content of this section has been taken directly from Schiller et al. [2] with minor linguistic 
changes. 
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radical reaction, a similar distribution of the final conversion degree could be 
observed: Highest in the pure resin layer and decreasing towards the fiber, with 
the lowest value directly at the fiber surface. Notwithstanding, it was found that 
the transition between the different layers was almost smooth. This is most 
likely due to the diffusion of some molecules, compensating for the change in 
the composition of sizing and resin layers. 


The introduction of the fiber surfaces reduced the diffusivity in the normal di- 
rection to the surface. This leads to a reduced overall conversion degree of the 
radical reaction of the final system, compared to a pure resin system. Moreover, 
this also results in locally varying conversion degrees and anisotropic radical 
polymerization at the fiber surface. A comparison of the results of this work to 
experiments would be highly desirable. Although various investigations of the 
fiber interface were conducted [68, 172, 173] a comparison of the results is not 
possible. These experiments are mainly focused on a mechanical characteriza- 
tion of the interphase based on micromechanical tests. This results in effective 
quantities for the whole fiber-sizing-resin interface. Furthermore, these exper- 
iments clearly show that a significant part of the failure of FRP is due to the 
interface between fiber and resin, including the sizing. Therefore, a better un- 
derstanding of the detailed processes during polymerization could also improve 
the design of such experiments. Nevertheless, any experiments investigating, 
for example, the diffusion of the sizing component during the reaction would 
allow a direct validation of the presented results. Moreover, this work now 
makes it feasible to include the mechanical testing of the final system within 
the MD. This is a non-trivial task, but it allows a direct comparison with the 
existing literature and could be the focus of a subsequent work. 


In other subsequent works, the complexity of the fiber sizing can be increased: 
Instead of only one coupling agent, several coupling agents or a different film 
former could be used. In addition, other additives of the sizing could be consid- 
ered. Also, the basic modeling of the condensation reaction could be extended. 
For example, a coarse-grained approach could provide detailed insight into the 
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behavior of the sizing during condensation, which could enhance the full atom- 
istic studies of the proposed approach. The assumption of the rigidity of the 
fiber surface could also be dropped, requiring the complex atomic structure 
of the fiber surface to be modeled. This could potentially improve the results, 
but would significantly increase the complexity. Last but not least, a further 
investigation ofthe generated system can be conducted. For example, an evalu- 
ation of the material properties of the system, e.g. the thermal or (visco-)elastic 
parameters, would be of great interest. In particular, an evaluation dependent 
on the distance to the fiber surface, as proposed in this work, would allow a 
deeper understanding of the sizing-resin interface. And eventually, complex 
properties such as a realistic interfacial fracture energy could be derived from 
such a system. 


Phase-field models for crack propagation in heterogeneous sys- 
tems? In this work, two different phase-field models for crack propagation 
in heterogeneous systems were introduced and compared: 


e A SCOP model, based on established approaches, which uses a single 
crack order parameter to account for damage. 


e A novel MCOP model, that introduces multiple order parameters in or- 
der to distribute the effective damage to the individual regions, modeled 
by its own set of order parameters. This results in multiple evolution 
equations, each of which has a constant crack surface energy. 


It was shown that the SCOP model is not able to reproduce the surface en- 
ergy of the sharp interface for an interfacial crack, especially when the same 
length parameters are chosen for the solid and crack problem. In comparison, 
the novel MCOP formulation avoids any errors in that case. Furthermore, the 
model demonstrated the same kinetics and crack profiles during the propaga- 
tion though an artificial interface, which confirms that multiple crack order 


?The content of this section has been taken directly from Schiller et al. [1] with minor linguistic 
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parameters can generate a continuously effective crack. An extension to an in- 
finitely tough interface indicated further problems with the SCOP model, as the 
crack is deflected from the interface before it hits the interface itself, and main- 
tains a certain distance from it. In contrast, the novel MCOP model demon- 
strates a more pronounced kink of the crack path at the interface and grows 
directly along the interface, for both sharp and diffuse interfaces. A quantita- 
tive comparison between the models for different angles of the interface and an 
analytical solution also showed a huge improvement in the modeling of crack 
propagation in heterogeneous materials: In particular, for high angles of the 
interface the SCOP model cannot replicate the analytical solution, where the 
MCOP model still shows very good agreement with it. 


The application of the novel model to FRP for unidirectional fibers in 2D and 
the extension to a 3D domain with isotropic fiber orientation distribution shows 
that the model is able to depict crack evolution in such a complex system, in- 
cluding crack nucleation and merging. Notwithstanding the fact that the crack 
phase-field model avoids established extensions, such as a tension-compression 
splitting or the removal of the interface parameter dependence, the many re- 
maining advantages of using multiple crack order parameters could be demon- 
strated. In future work, a tension-compression split, for example based on the 
work of Storm et al. [148], could further improve the model. In addition, a 
combination with a solid-solid phasefield transition model and plasticity could 
allow the study of other material systems, e.g., crack evolution during marten- 
sitic phase transformation [174]. The inclusion of an interfacial crack resis- 
tance could also allow for a more sophisticated simulation of FRP failure, by 
enabling realistic fiber debonding. 


Considering mechanical jump conditions in the MCOP model? In 
this work a novel MCOP model is proposed for fracture in heterogeneous ma- 
terials. Despite the improvement in the qualitative and quantitative prediction 


3This section is based on the work of Schiller et al. [3]. Minor linguistic changes and additions 
have been made. 
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of crack paths for such systems, the model has some limitations. Therefore, an 
extension of this model was proposed in this work and investigated: Instead 
of a basic Voigt-Taylor homogenization scheme, the approach is extended to 
consider mechanical jump conditions based on Schneider et al. [61]. First, the 
single inclusion problem was introduced to demonstrate the difference in elas- 
tic driving forces of the two schemes. In addition, the final crack paths due to 
eigenstrains mimicking chemical shrinkage were discussed for the single inclu- 
sion problem for both schemes as well as the single-well and single-obstacle 
schemes. In this context, the jump condition scheme resulted in a better agree- 
ment with the expected fracture behaviour. 


An exemplary FRP system was introduced to investigate the behavior of both 
schemes. Therefore, the crack resistance of the glass fiber was varied. The 
Voigt-Taylor homogenization scheme failed to predict reasonable crack paths 
for contrary contrasts of elastic modulus and crack resistance. Instead, the 
fiber failed as well, resulting in different paths. In contrast, when mechanical 
jump conditions are considered, the model yields the same final crack path for 
all crack resistances presented, since the mechanical driving force for crack 
propagation is modeled more independently of the elastic properties of other 
physical domains. Moreover, this behavior could also observe in the stress- 
strain curves. While the Voigt-Taylor scheme fails at different loading points, 
the novel scheme shows a negligible scatter. 


Based on this work, more extensive simulation studies on the failure mecha- 
nism of FRP materials, as well as other complex systems such as polycrys- 
talline materials, e.g., in hydrothermal environments [175], solid oxide fuel 
cells among others, can be conducted. Combined with experimentation, this 
offers the opportunity to improve a broad range of engineering applications. 


Crack propagation during curing in fiber-reinforced thermosets 
In order to predict the crack formation of FRTS during curing, a method to 
generate virtual fiber-reinforced volume elements was introduced in Chapter 4. 
For this purpose, an approach based on MD was used to consider not only 
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straight fibers, but also long curved fibers. In the following Chapter 5, the 
curing process, which was previously modeled on the atomic length scale, was 
simulated on the micro length scale using continuum mechanics. Incorporat- 
ing thermal and chemical strains, the resulting residual stresses were discussed 
for different levels of abstraction of the microstructure of FRTS. Finally, in 
Chapter 7, the various aspects of this work were brought together. Based on 
the novel crack propagation model for heterogeneous systems and the curing 
model for FRTS, crack formation during manufacturing was simulated. De- 
spite the lack of possibility to compare the results with experiments and the 
difficulty to assess the quantitative results, the cracking mechanisms of this 
material class could be observed. It was demonstrated that the presented ap- 
proach is able to show crack formation during curing of glass fiber-reinforced 
UPPH resin. For further work, such a model could benefit from incorporating 
the viscoelastic material behavior of the resin, despite the increase in compu- 
tational cost and lack of material properties during the various stages of the 
curing process. Based on further experimental investigation and incorpora- 
tion of crack nucleation and interfacial crack propagation, a more quantitative 
prediction of crack formation during the curing process seems feasible and 
would be highly desirable for understanding this material class and its failure 


mechanism. 
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A.1 Intermolecular interaction potentials 


The intermolecular interaction potential energy function Ui", cf. (2.64), con- 
sists of several contributions introduced below and are based on the COMPASS 
force field [101-103]. Note that the individual parameters in these contribu- 
tions are usually dependent on the atom types, but for the sake of clarity addi- 
tional sub or superscripts are omitted. 


Bonds The bonds contribution consist of a potential with harmonic and an- 
harmomic terms 


UR =K} (rij— ro) +K? (rij— ro)? +K? (rij — ro) (A.1) 


with the parameters K}, K}, K} and the equilibrium bond distance ro. 


Angles The angle contribution can further decomposed in 


with a angle potential, again with harmonic and anharmomic terms 


4. = K3 (Oije — 00)? +K (Oije — 80)” + K3 (0; — 8) (A.3) 
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with the parameters K},K3,K5 and the equilibrium angle &%. The bond-bond 


term follows by 
Up = om? (nj _ rt) (r; k~ ar (A.4) 


with the parameters M®? and the distances r; and ae In addition, bond-angle 
term 


Ups = NP (nj = ri") (6; jx — 00) +N2 (rx — i") (6; jx — 90) (A.5) 


with the parameters N ba, ypa and the distances pa and rei 


Dihedrals The dihedrals contribution can further decomposed in 


UP 


l 


Da = Uf + Uie +U pia + Ui + Ui + Up > (A.6) 


with a dihedral potential 
d ee d 
Usa = È Ki (1— cos (nd; - 08) ) (AN) 
n=1 


with the parameters K d KÌ, KÌ and angles of, os, os. The middle-bond-torsion 
term follows by 


Um = (rx — ry) (Aı cos (Qj jx) +Ao2cos (26; jx) +A3cos (36; jx) ) ; 
(A.8) 
with the parameters A;,A2,A3 and the distance pe Moreover the end-bond- 
torsion term 


U = (nj — ii") (Bı cos (Qij) + B2 cos (26; jx) + B3 cos (36; x) ) + (A.9) 


(mr? ye cos (Pijer) +C2cos (20; ja) + C3 cos (3;;x1)) , (A.10) 
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with the parameters B1, B2,B3,C1,C2,C3 and the distance r‘®',r$°' and the 
angle-angle-torsion term 


ikl = (jx _ 8°) (Dı cos (Oj 1) + D2 cos (26; jx) + D3 cos (30; jx) ) + 
(A.11) 


(Aji — 05°) (E1 cos (Qij) + E2 cos (2; x1) + E3 cos (39; jx1) ) ; 
(A.12) 


with the parameters D1, D2,D3, E1, E2, E3 and the angles 6}, 05°'. In addition, 
the angle-torsion follows by 


Uli = M* (Oije — 01") (Oj — 5") cos (Qijki), (A.13) 


with the parameters M“ and the angles 67", 03" and the bond-bond-13 term 
Bee — Npb13 (nj = P’) (ru = ae (A.14) 


Nb! ae deed 


with the parameters and the distances r? 


Impropers The impropers contribution can further decomposed in 


Un = Uh igs + Un (A.15) 


with the improper term 


Al 2 
Uiju = K' G (Xijki + Xkjti + Kıjik) — x) (A.16) 
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with the parameters Kİ and the equilibrium angle 7%. The angle-angle term 
follows by 


By =M¥ (Oije — 01°) (Ok; — 05°) + (A.17) 
M3" (Oije — 03°) (Oii — 03°) + (A.18) 
M3" (O;jı — 03°) (Okjı — 03°) (A.19) 


with the parameters Mj", M5*, M3" and the angles 07°, 05°, 03". 


A.2 Two-stage polymerization! 


Benzene distribution In order to be able to investigate the spatial distri- 
bution of the various constituents, a reasonable approach must be chosen. Es- 
pecially molecules like (P-)MDI or the unsaturated polyester consist of numer- 
ous atoms, which can be distributed widely over the system. Reducing them to 
their center of mass would not result in a useful distribution. Instead, specific 
functional groups are chosen. For example, all molecules (except the coupling 
agent) contain benzene rings. In addition, they can be distinguished by their 
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Figure A.1: Spatial distribution of the benzene of styrene during radical polymerization. 


'The content of this section has been taken directly from Schiller et al. [2] with minor linguistic 
changes. 
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Figure A.2: Spatial distribution of benzene of (P)-MDI during polyurethane reaction (left) and 
radical polymerization (right) for different conversion degrees. 
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Figure A.3: Spatial distribution of benzene of the unsaturated polyester during polyurethane reac- 
tion (left) and radical polymerization (right) for different conversion degrees. 


uniquely bonded atoms. This allows easy tracking of all benzene rings across 
the system and during polymerization reactions. 


Figure A.1 shows the styrene benzene during radical polymerization. During 
the reaction, the small styrene molecules level out the initial uneven distribu- 
tion due to diffusion. The benzene distribution of the larger molecules such 
as (P-)MDI and the unsaturated polyester is shown in Figure A.2 and Fig- 
ure A.3, respectively. Despite the fact that they also exhibit some additional 
artificial charges during the polyurethane reaction, significant diffusion is not 
observed. In contrast, they show diffusion during the radical polymerization, 
as they smooth the initial uneven distribution. 
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Figure A.4: Spatial distribution of silicon atoms of y-MPS molecules during polyurethane reaction 
(left) and radical polymerization (right) for different conversion degrees. 


Silicon distribution of coupling agent For the distribution of the cou- 
pling agent, hence y-MPS, the silicon atom of the molecules is used to rep- 
resent the position of the free hydrolyzed as well as the coupling agent pre- 
attached to the fiber surface. Figure A.4 shows this distribution during the 
polymerization reactions. As before, no significant diffusion can be observed 
during the polyurethane reaction. Nevertheless, the pre-attached layer with the 
very high peak at the fiber surface as well as the absence of Y-MPS in the pure 
resin layer can be observed. During the free radical polymerization, as with the 
other constituents, diffusion takes place, smoothing out the initial distribution. 


A.3 Phase-field modeling of crack 
propagation? 


Staggered iterations study For a numerical investigation of the number 
of staggered iterations during a crack propagation simulation, the setup of an 
artificial sloped interface, cf. Section 6.4.2, is chosen. For this purpose, the 
number of iterations is depicted in Figure A.5. The cumulative iterations over 


?The content of this section has been taken directly from Schiller et al. [1] with minor linguistic 
changes. 
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Figure A.5: Behavior of the staggered scheme for the artificial sloped interface: (a) Cumulative 
iterations during the crack growth for the SCOP and MCOP models. (b) Amount of 
inner iterations for an exemplary discrete crack propagation step. 


the crack growth are shown in Figure A.5a. For this purpose, all inner staggered 
iterations are summed up until a discrete crack propagation occurs. Thereby, 
the SCOP and MCOP models show similar behavior. The number of iterations 
is approximately equal. Only near the solid interface the MCOP requires more 
iterations. This is most likely due to the implicit transition of the crack from 
one order parameter to the other. The overall trend towards fewer iterations, 
like the change in slope, is likely due to the way the stress boundary condtion 
is imposed. In addition, Figure A.5b plots the number of staggered iterations 
against temporal iterations for a discrete crack propagation step. Also here, 
SCOP and MCOP show the similar behavior: After a discrete crack propaga- 
tion, the first iterations consist of an increased number of staggered iterations, 
while subsequently, the iterations are lower, but increase slightly again to the 


next discrete crack propagation. 
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Figure A.6: Mesh convergence study of the SCOP model for the artificial sloped interface: The 
mesh sizes at the crack tip Ax are varied. In addition, the background was also refined 
and the adaptive mesh refinement parameters were chosen more cautiously (+). 


Mesh convergence study For an investigation of mesh convergence for 
the artificial sloped interface example (cf. Section 6.4.2), the SCOP model and 
the infinitely hard upper region are chosen. As in this case the highest gradients 
occur in the solution fields and can therefore be assumed to be the most chal- 
lenging problem to discretize. Figure A.6 shows the results of the SCOP model 
of Section 6.4.2. In addition, the mesh size Ax at the crack tip is varied, where 
Axo is the size in of the previous results. Furthermore, the underlying coarse 
mesh is significantly refined, and the adaptive mesh refinement parameters are 
changed to increase the area where Ax applies, denoted by + in Figure A.6. 
Thereby, none of the results yields large variance, therefore the discretization 
chosen in this work can be assumed to be representative. 
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